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In this work we present 3D numerical relativity simulations of thick accretion disks around tilted 
Kerr black holes. We investigate the evolution of three different initial disk models with a range of 
initial black hole spin magnitudes and tilt angles. For all the disk-to-black hole mass ratios considered 
(0.044 — 0.16) we observe significant black hole precession and nutation during the evolution. This 
indicates that for such mass ratios, neglecting the self-gravity of the disks by evolving them in a 
hxed background black hole spacetime is not justified. We hnd that the two more massive models 
are unstable against the Papaloizou-Pringle (PP) instability and that those PP-unstable models 
remain unstable for all initial spins and tilt angles considered, showing that the development of the 
instability is a very robust feature of such PP-unstable disks. Our lightest model, which is the most 
astrophysically favorable outcome of mergers of binary compact objects, is stable. The tilt between 
the black hole spin and the disk is strongly modulated during the growth of the PP instability, 
causing a partial global realignment of black hole spin and disk angular momentum in the most 
massive model with constant specific angular momentum 1. For the model with non-constant l- 
profile we observe a long-lived m = 1 non-axisymmetric structure which shows strong oscillations of 
the tilt angle in the inner regions of the disk. This effect might be connected to the development of 
Kozai-Lidov oscillations. Our simulations also confirm earlier findings that the development of the 
PP instability causes the long-term emission of large amplitude gravitational waves, predominantly 
for the I — m — 2 multipole mode. The imprint of the BH precession on the gravitational waves from 
tilted BH-torus systems remains an interesting open issue that would require signihcantly longer 
simulations than those presented in this work. 

BAGS numbers: 04.25.D-, 04.30.Db, 95.30.Lz, 95.30.Sf, 97.60.Lf 97.10.Gz 


I. INTRODUCTION 

Stellar mass black hole-torus systems are believed to 
be the end states of binary neutron star (NS-NS) or black 
hole neutron star (BH-NS) mergers, as well as of the 
rotational gravitational collapse of massive stars. BH- 
torus systems emit gravitational waves (GWs), which 
may eventually provide the direct means to study their 
actual formation and evolution and prove the current hy¬ 
potheses that associate them to gamma-ray burst (GRB) 
engines mm- Such a possibility is out of reach to electro¬ 
magnetic observations due to their intrinsic high density 
and temperature. For an overview of the event rate esti¬ 
mates of NS-NS and BH-NS mergers that are observable 
with initial and advanced LIGO see e.g. [SE]. 

Our theoretical understanding of the formation of BH- 
torus systems and their evolution relies strongly on nu¬ 
merical work. In recent years a significant number of nu¬ 
merical relativity simulations have shown the feasibility 
of the formation of such systems from generic initial data 
(see e.g. [GHIO] for recent progress). In particular, the 
3D simulations of [6] (see also references therein) have 
shown that unequal-mass NS-NS mergers lead to the self- 
consistent formation of massive accretion tori (or thick 
disks) around spinning black holes, thus meeting the nec¬ 
essary requirements of the GRB’s central engine hypoth¬ 


esis. However, if the energy released in a (short) GRB 
comes from the accretion torus, the BH-torus system has 
to survive for up to a few seconds m- Any instability 
which might disrupt the system on shorter timescales, 
such as the runaway instability m or the Papaloizou- 
Pringle instability (PPI hereafter) [13], could pose a se¬ 
vere problem for the prevailing GRB models. 

Using perturbation theory, Papaloizou and Pringle 
found that tori with constant specific angular momentum 
are unstable to non-axisymmetric global modes. Such 
modes have a co-rotation radius within the torus, located 
in a narrow region where waves cannot propagate. Waves 
can still tunnel through the co-rotation zone and interact 
with waves in the other region. The transmitted modes 
are amplified due to the feedback mechanism provided 
by the reflecting boundaries at the inner and outer edges 
of the torus. The manifestation of the PPI, as [T4| first 
showed numerically through hydro dynamical simulations 
in the fixed background metric of a Schwarzschild BH, 
is in the form of counter-rotating epicyclic vortices, or 
“planets’”, with m planets emerging from the growth of 
a mode of order m. Moreover, BH-torus systems are char¬ 
acterized by the presence of a cusp-like inner edge where 
mass transfer driven by the radial pressure gradient is 
possible. If due to accretion the cusp moves deeper inside 
the disk material, the mass transfer speeds up leading to 
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a runaway process that destroys the disk on a dynamical 
timescale (see e.g. m for test-fluid simulations in gen¬ 
eral relativity of the occurrence of this instability and [16] 
for axisymmetric simulations where self-gravity was first 
taken into account). In most recent numerical relativ¬ 
ity simulations [aiaiiniin] the BH-torus systems under 
consideration did not manifest signs of dynamical insta¬ 
bilities on short dynamical timescales, because the non¬ 
constant angular momentum profiles of the massive disks 
that formed in the simulations seem to make them stable 
against the development of the runaway instability. How¬ 
ever, on longer timescales non-axisymmetric instabilities 
can easily develop [13 HH. 

Korobkin et al m have studied non-axisymmetric 
instabilities in self-gravitating disks around BHs us¬ 
ing three-dimensional hydro dynamical simulations in full 
general relativity. Their models incorporate both mod¬ 
erately slender and slender disks with disk-to-BH mass 
ratios ranging from 0.11 to 0.24. While no sign of the run¬ 
away instability was observed for these particular models, 
unstable non-axisymmetric modes were indeed present. 
More recently [20| observed that by a suitable choice 
of model parameters, namely constant angular momen¬ 
tum tori exactly filling or slightly overflowing their Roche 
lobe, a rapid mass accretion episode with the character¬ 
istics of a runaway instability sets in. The astrophysi- 
cal significance of such fine-tuned models is uncertain as, 
e.g. they do not seem to be favored as the end-product of 
self-consistent numerical relativity simulations of binary 
neutron star mergers 1313113 HI]- Correspondingly, the 
3D general relativistic numerical simulations of BH-torus 
systems of [T9| showed that an m = 1 non-axisymmetric 
instability grows for a wide range of self-gravitating tori 
orbiting BHs. Their models included torus-to-BH mass 
ratios in the range 0.06 — 0.2 and both constant and non¬ 
constant angular momentum profiles in the disks. [T9| 
found that the non-axisymmetric structure persists for a 
timescale much longer than the dynamical one, becoming 
a strong emitter of large amplitude, quasi-periodic GWs, 
observable by forthcoming detectors. 

The vast majority of numerical relativity simulations 
to date that have studied the formation of BH-torus sys¬ 
tems have led to systems in which the accretion torus is 
aligned with the equatorial plane of the central Kerr BH. 
There are, however, several ways of arriving at tilted ac¬ 
cretion disks around Kerr BH. Possible scenarios include 
asymmetric supernova explosions in binary systems [22| 
or via BH-NS mergers in which the BH spin is misaligned 
with the orbital plane of the binary. See also [23l |24| for 
arguments that most BH-torus systems should be tilted. 
Recently, full GRHD simulations of BH-NS mergers with 
tilted initial BH spins have been performed by [254127] . 
which can result in tilted accretion tori around the BH 
remnant. In these simulations massive remnant disks only 
formed for high initial BH spins. As shown in [28], the 
disk mass in BH-NS mergers increases with increasing 
initial BH spin and decreases with increasing initial BH 
mass. This is due to the size of the ISCO (innermost cir¬ 


cular stable orbit) of the BH in the merger. The ISCO 
grows with BH mass and decreases with the spin magni¬ 
tude of the BH. If the ISCO is large enough, the NS will 
be ’swallowed’ entirely by the BH before being tidally 
disrupted, leaving no accretion torus behind. In order to 
estimate disk masses of resulting from BH-NS mergers, 
one needs thus an estimate for the initial BH masses in 
these system. One such method, via population synthe¬ 
sis considerations, favor larger BH masses [29413T] . with 
a peak around 8 Mq and a mass gap between the light¬ 
est expected BHs and NS masses. This means that these 
massive BH would need very large initial spins in order 
to be able form massive remnant disks after the BH-NS 
merger. 

General relativity hydrodynamics and MHD simula¬ 
tions of such kind of tilted thick accretion disks around 
BHs have been performed by Fragile and collaborators 
within the test-fluid approximation (see e.g. [32l l33]b 
These authors have carried out an extensive program to 
study both the dynamics and observational signatures of 
thick accretion tori around tilted Kerr black holes. 


Motivated by these previous works we present in this 
paper an extended set of numerical relativity simulations 
of tilted accretion disks around Kerr BHs in which the 
self-gravity of the disks is taken into account. In this first 
study of the effects of the disk self-gravity in tilted BH- 
torus systems, we have chosen low mass central black 
holes, and small spins that are unlikely to produce mas¬ 
sive tilted disks in astrophysical mergers for the reasons 
outlined above. As a first study, we want to investigate 
the effects of the disk self-gravity around tilted BH with 
our existing initial data, where some of our initial mod¬ 
els were known to be unstable against the PPI in the 
untilted case [18]. We plan on studying these tilted BH- 
torus systems with more astrophysically realistic data in 
future works. Our simulations incorporate the effects of 
the self-consistent evolution of BH-torus systems with 
misaligned spins via the solution of the full 3-1-1 Ein¬ 
stein equations coupled to the hydrodynamics equations. 
Due to the Lense-Thirring effect, the inner regions of 
the disks start precessing and might become twisted and 
warped affecting their dynamical behavior (for a defini¬ 
tion of twist and warp see e.g. [34]). Those inner regions 
might also be forced into the equatorial plane of the BH 
due to the so-called Bardeen-Petterson effect [35]. Our 
simulations allow us to investigate these issues and find 
out how the dynamics of the tilted disks affect the cen¬ 
tral BH, assessing in turn the importance of neglecting 
the disk’s self-gravity on the overall dynamics of these 
systems. 


This paper is organized as follows: In Section [TT| we 
summarize the mathematical and numerical framework 
we employ in our numerical study, explaining in detail 
the methods we use for the spacetime and matter evo¬ 
lution. We also describe several diagnostics we use to 
monitor and analyze the simulations. In Section HI we 


describe the setup of our initial data, both in the un¬ 
tilted and tilted cases. We also provide a table which 
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contains the key parameters of the models considered in 
this study. Our results are presented and discussed in de¬ 
tail in Section [TV] Finally, we present our conclusions in 
Section [V| Throughout the paper, we use units in which 
c = G = Mq = 1 where c, G and Mq are the speed 
of light, the gravitational constant and the solar mass, 
respectively. Latin indices run from 1 to 3 while Greek 
indices run from 0 to 3. 


II. MATHEMATICAL AND NUMERICAL 
FRAMEWORK 

A. Spacetime evolution 

The simulations reported in this paper are per¬ 
formed using the publicly available Einstein Toolkit (ET) 
code [Ml EZl- The ET is a code for relativistic as¬ 
trophysics simulations, which uses the modular Cactus 
framework [38] (consisting of general modules commonly 
called Thorns’) and provides adaptive mesh refinement 
(AMR) via the Carpet driver [39]. The set of com¬ 
plex time-dependent partial differential equations is in¬ 
tegrated using the met ho d-of-lines where the time inte¬ 
gration is performed using a 4th order Runge-Kutta al¬ 
gorithm. The toolkit provides numerical solvers and the 
necessary infrastructure to evolve the Einstein equations, 

=8 7rT^^, (1) 

coupling the evolution of spacetime to the dynamic of 
the matter content through the stress-energy tensor 
The methods to evolve the Einstein equations are based 
on the 3+1 ADM |40] splitting of spacetime, in which the 
metric becomes 

ds^ = (—+ 2/3^ dt dx'^ + dx'^ dx^ , (2) 

where a is the lapse function, the shift vector, and 
^ij is the spatial metric tensor. The extrinsic curvature 
tensor Kij is constructed from the spatial metric as 
follows: 

Kij = ~ ^ 13 ) Jij 5 (3) 

where Cj^ stands for the Lie derivative with respect to 
the shift vector. 

The left hand side of the Einstein field equations is 
evolved using the McLachlan code miiii], which solves 
a conformal-traceless “3 + 1” formulation of the Einstein 
equations known as BSSN [434145] . 

The BSSN evolution equations are evolved using a 
fourth-order accurate, centered finite-differencing oper¬ 
ator, while the advection terms for the shift vector 
are evolved with a fourth-order upwind stencil. We ap¬ 
ply fifth-order Kreiss-Oliger Dissipation to all space- 
time variables to achieve overall fourth-order for the 
spacetime evolution. Using the BSSN evolution system 
and the gauge conditions described in Appendix the 


McLachlan code can evolve stably a single puncture for 
several light crossing times, see e.g. m 

We apply Sommerfeld outgoing boundary conditions 
for all the BSSN evolution variables with an extrapola¬ 
tion to include the part of the boundary conditions that 
does not behave like a simple outgoing wave [46] . 

B. Matter Evolution 

The evolution of the hydrodynamics is performed by 
GRHydro [37l iTfl l48]. that solves the general relativis¬ 
tic hydrodynamics equations in flux-conservative form, 
in the so-called Valencia formulation |494[5T] , using high- 
resolution shock-capturing (HRSC) methods. The GRHD 
evolutions equations are described in Appendix [B| 

As a grid-based code GRHydro cannot evolve regions 
without matter content due to the breakdown of the EOS 
and evolution equations in regions of zero rest-mass den¬ 
sity. As is customary, a low density atmosphere (typically 
several orders of magnitude below the maximum density 
of the initial data) fills those grid points that should be¬ 
long to vacuum regions. In all our simulations we set 
the density of the atmosphere to 10“^ times the initial 
maximum rest-mass density in the torus. Although fully 
evolved, the dynamical effects of the atmosphere regions 
can be neglected. 

We use the piecewise parabolic method (PPM) [52] to 
reconstruct the primitive variables at the cell interfaces 
and Marquina’s flux formula [53] [54] to compute the nu¬ 
merical fluxes. Differently from the original implementa¬ 
tion reported in Ref. [Ml, we reconstruct the quantities 
Wv'^ instead of the three-velocities v'^. This guarantees 
that the velocities reconstructed at the cell boundaries 
remain subluminal even under extreme conditions like 
those encountered near the apparent horizon (AH) of the 
BH. 


C. Diagnostics 

We measure the BH mass during the evolution using 
the AHFinderDirect thorn which implements the AH 
finder described in [56]. The BH spin magnitude and 
direction are measured using the QuasiLocalMeasures 
thorn EH EH] (see also the review of quasi-local meth¬ 
ods in (Ml). We measure the spin direction using the so- 
called flat-space rotational Killing vector method of m, 
which can be derived using Weinberg’s pseudotensor in 
Gaussian coordinates and is equal to the Komar angu¬ 
lar momentum integral when the latter is expressed in a 
foliation adapted to the axisymmetry of spacetime m- 
To track the moving location of the puncture during the 
evolution, we use the PunctureTracker thorn. Gravita¬ 
tional waves are extracted using a multipole expansion of 
the Weyl scalar at different radii. The thorns Multipole 
and WeylScal4 compute the respective quantities. We 
calculate the mass accretion rate M as the instantaneous 
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FIG. 1: The figure visualizes the various vectors and projec¬ 
tions we use to calculate the twist and tilt in the disk. Note 
that all vectors in the hgure have been scaled to unity for 
better visualization. 


flow of matter through the AH, using the Outflow thorn. 
To quantify the growth of non-axisymmetric modes in the 
disk, we use the thorn GRHydroAnalysis which computes 
3D Fourier integrals of the density, described in [62] . The 
amplitude of the m-th mode is given by: 


Dr, 


= A. 


(4) 


Since the computational domain is very large (see sec¬ 
tion III C| below), the total mass of the atmosphere can 
be a non-negligible fraction of the total rest-mass of the 
torus. We therefore have to ignore cells corresponding 
to the atmosphere when integrating the torus rest-mass 
density and in the calculations of the non-axisymmetric 
modes in the disk. Despite the non-negligible mass of the 
atmosphere, this mass is not dynamically important be¬ 
cause it is spatially homogeneous (and therefore its grav¬ 
itational pull on the central BH-torus system is zero) and 
its coordinate velocity is set to zero. 


D. Twist and tilt angles 


vector, we define two angles, the twist a{r) 

(j{r) = Z(Sbh X Sxy-90, P(JDisk(^), Sbh)) , (5) 

and the tilt v{r) 

v{r) = A(Sbh, JDisk(^)) , (6) 

where 

P(a,b) = a-(7) 

is the projection of vector a onto the plane with normal 

b. 

The vector Sbh x Sxy -90 is constructed in the following 
way: We project the BH spin Sbh onto the x^-plane, 
Sxy = P(Sbh, z) and then rotate the resulting vector Sxy 
by 7r/2 about the 2 :-axis. The cross product of Sbh and 
Sxy -90 then lies in the equatorial plane of the BH and 
the twist cr(r) at t = 0 is 0 throughout the disk. A sketch 
showing the construction of these vectors and twist and 
tilt angles is provided in Fig. 

The twist <j(r) is then a measure of how much the 
angular momentum vector of each annulus has precessed 
in the (dynamically changing) equatorial plane of the BH, 
while the tilt z^(r) gives the angle between the BH spin 
and the angular momentum vectors for each annulus. The 
disk is said to be twisted if the twist becomes a function 
of radius, a = cr(r), and said to be warped if the tilt 
becomes a function of radius, u = u{r). 

The angular momentum vector for each annulus is cal¬ 
culated following the procedure outlined in the section 
on spin in [63] (pages 46f.). We calculate JDisk(^) from 
the anti-symmetric angular momentum tensor 

L^^•' = j , ( 8 ) 

in the following way: 

^ j31 ^ ]^12 _ (g) 

The relevant components of the stress-energy tensor 
are 


= ph 




+ pT 

~ -I- r, 5 


( 10 ) 


In order to keep track of the response of the tilted ac¬ 
cretion disk and to check for Lense-Thirring precession 
and the occurrence of the Bardeen-Patterson effect, we 
measure two angles, the twist (precession) and tilt (incli¬ 
nation) of the disk [34| . We closely follow the description 
of the angles given in [32], with a modification in the 
calculation of the twist, and furthermore adapt the way 
they are calculated. 

The main idea is to split the disk up into a series of 
annuli and calculate the angular momentum vector of the 
matter, JDisk(^), for each individual annulus. Using this 


where Latin indices run from 1 to 3. 

We note that our procedure differs slightly from that 
given in [32], where the components of the disk shell an¬ 
gular momentum vectors were computed using the intrin¬ 
sic spin as follows: 

Sa = l , ( 11 ) 

where e^^^ys is the 4-dimensional Levi-Civita symbol and 
W = p^i-ppp^)^ is the total 4-velocity of the system, 
with being the 4-momentum. 
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Using the extrinsic spin, the contribution to the total 
angular momentum stemming from choosing a reference 
point is accounted for in the linear momentum of the 
center of mass of the system. We choose to use the angu¬ 
lar momentum tensor calculated about the center 

of the BH instead, using formula (§ in the calculation of 
cr(r) and u{r). We therefore pick a reference frame cen¬ 
tered on the BH, as it is moving around the grid and we 
are interested in the total angular momentum of the disk 
shell about the BH center. We note that the calculation 
of the twist and tilt angles is in fact gauge dependent. 
We have checked the impact of the gauge choice on the 
twist and tilt angles by measuring the angles for the ini¬ 
tial data with different lapse and shift profiles to find 
that the angles are calculated correctly (as the initial 
tilt and twist profiles are parameters of our initial data), 
which also serves as an important consistency check for 
the correct measure of the initial angles. In the actual nu¬ 
merical implementation, we do not use the inverse cosine 
to calculate the angles, as the function becomes very in¬ 
accurate when the vectors become close to being parallel 
or anti-parallel, instead we use the atan2 function which 
is free of this deficiency. The tilt u{r) is then calculated 
using the following formula: 


v{r) = tan ^ |Sbh x JDisk(^)| , Sbh • J 




( 


where the hat indicates unit vectors. 

For the twist cr(r), we need to be more careful. In the 
formula above, we implicitly assume that v{r) < tt , V t, 
because Eq. (HI will always return angles < tt. Although 
it does not generally make much sense to consider an¬ 
gles > TT between two 3D vectors, because there is an 
up and down direction, we are nevertheless interested in 
directional (signed) angles > tt for the twist. The rea¬ 
son is that we want to be able to track cumulative twists 
larger than tt. In order to calculate the directional angle 
cr(r) from the reference vector Sbh x Sxy -90 to the tar¬ 
get vector P(JDisk(^), Sbh) in a fixed sense of rotation, 
we need a third reference vector that always lies above 
the plane spanned by the two original vectors. As the 
plane in which both vectors live is constructed to be the 
equatorial plane of the BH, we can therefore choose this 
vector to be Sbh- This information allows us (because 
the cross product contains this information in the sign 
it returns) to calculate twist angles in the range [0, 27r] 
with the following formula: 

(j{r) = tan“^ [ |Jbh x JDisk(^)| , Jbh Joisk (r)] (13) 

and then select the directional angle by: 

(rr.\ — f • ^BH • (Sbh X JDisk(^)) > 0 

^ \ 27r — cr(r) : Sbh * (Sbh x JDisk(^)) < 0 

We also compute the instantaneous precession and nu¬ 
tation rates for both the BH spin vector S and the total 
disk angular momentum vector Joisk (which we calculate 
by summing the components of all shells). For the BH, 
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FIG. 2: Innermost region of the initial rest-mass density pro¬ 
file along the x-axis for model ClBaObO. The disk ends at 3.03 
and the AH is initially located at x = 0.46. 


we compute its precession and nutation about the z-axis, 
while the total precession and nutation of the disk angu¬ 
lar momentum vector Joisk is calculated about the BH 
spin axis. 


III. INITIAL DATA AND SETUP 


Our initial setup is a thick, self-gravitating axisym- 
metric accretion disk in equilibrium around a rotating 
BH. Such a system is built following the approach laid 
out in |64] that we briefly describe below. The solver to 
build the initial data (ID) first computes models of self- 
gravitating, massive tori around non-rotating BHs. Then, 
for our simulations of tilted disks around rotating BHs 
we retain the hydrodynamical content of a model and re¬ 
place the spacetime by a tilted Kerr spacetime in quasi¬ 
isotropic coordinates. The generated initial data is then 
interpolated onto Cartesian coordinates and evolved with 
the Einstein Toolkit. We are not aware of a method 
to generate self-consistent data for self-gravitating tilted 
disks around Kerr black holes and therefore resort to 
this method. In order to keep the constraint violations 
as low as possible, we use disk models with rest masses 
of only a few percent of the central BH mass and with 
small to moderate BH spins. When replacing the orig¬ 
inal spacetime of the initial data with the tilted Kerr 
spacetime in the computational domain, both the torus 
rest-mass Mq and the torus gravitational mass Mr will 
change somewhat, due to the fact that the volume ele¬ 
ment of the spacetime where g is the determinant 

of the four-dimensional metric g^iy, as well as the lapse, 
have changed. In order to assess the difference in Mq and 
Mt due to the replacement of the spacetime we calculate 
both quantities as in [64], as follows: 

Mq = f pu^y/^d^x= f pd^x, (14) 
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and 

Mt = I (-2 T! + T;) a (15) 

We find the largest difference in Mq and Mr in model 
ClBa03b30, where we obtain the following fractional dif¬ 
ferences: A Mo ~ —3% and A Mt ~ —3%. 


A. Self-gravitating accretion disks 


Our ID are self-gravitating disks around a 
Schwarzschild BH in quasi-isotropic (QI) coordi¬ 
nates [64]. The ID is fully described by the four metric 
potentials A(r, ^), a{r,0) and uj{r,0) that 

characterize the metric of a stationary, axisymmetric 
spacetime in spherical polar QI coordinates, 

ds^ = dt^ {df‘^d<p^)-\-^^f‘^sin^0{d(l)—ujdt)‘^ , 

(16) 

where the isotropic radius r is given in terms of the areal 
radius by 

f = i — M + \/r‘^ — 2Mr^ , (17) 


where M is the mass of the BH. The ID also provides 
the pressure p and orbital angular velocity O = /vd ^ 
measured by an observer at infinity at rest. The evolution 
time is measured in terms of the initial orbital period 
at the radius of the initial maximum of the rest-mass 
density, pc- These two quantities, together with the EOS, 
are sufficient to obtain the remaining hydro dynamical 
quantities. The ID is c onst ructed using the polytropic 
EOS described in Eq. (B9). The components of the 3- 
velocity in Cartesian coordinates are given by 


-O 


a 


-n 


a 


,0 


(18) 


where a is the lapse function. The Lorentz factor is cal¬ 
culated directly from the 3-velocity 

W = , ^ ^ . (19) 

^l-(w-n)2^(x2+y2) 


The matter content of the ID only fills the computa¬ 
tional domain up to the event horizon of the central BH. 
In our simulations we do not excise the BH but rather 
treat it as a puncture. We noticed that when evolving the 
BH as a puncture with matter content without excision, 
the simulations were not long-term stable and failed at 
early stages of the evolution due to errors in the mat¬ 
ter evolution at the location of the puncture. The rea¬ 
son is a very rapid pile up of matter at this location, 
which leads to non-physical values of the hydrodynami- 
cal variables at the puncture. These lead to failures in the 
GRHydro evolution code during the conversion from con¬ 
served to primitive variables. In [65] a successful method 


to evolve hydrodynamics in the presence of punctures 
was presented. This method checks for unphysical values 
of the hydro dynamical variables at the immediate vicin¬ 
ity of the puncture and resets them to physical values. 
In practice it checks for the positivity of the conserved 
energy (r > 0) and an upper bound for |5'p: 

Y^SiSj < t{t -1- 2p). (20) 

In our simulations, we use a similar treatment inside the 
BH horizon, namely we specify a fraction of the volume of 
the AH in which we reset the matter fields to atmosphere 
values and the stress-energy tensor to zero, which is 
essentially the same technique as described in [66] . In this 
way, we avoid using a moving excision boundary for the 
hydrodynamical variables as in M- Erom the point of 
view of the evolution it is safe to do this, as the evolution 
of the hydrodynamics inside the horizon cannot influence 
the evolution outside the horizon. In fact, we checked that 
the fraction of the AH in which we apply such atmosphere 
resetting has no effect on the evolution, as expected. We 
checked this by applying the atmosphere reset for differ¬ 
ent fractions of the AH. In all cases, the evolution out¬ 
side the AH was unaffected, confirming that the region 
inside the AH is causally disconnected. This has been 
employed for instance in the so-called “turduckening “ of 
initial BH data mi- Eor the practical implementation of 
this approach we use a spherical surface that contains 
the shape of the AH and apply the atmosphere method 
in a fraction of the minimum radius of that surface. Our 
procedure is illustrated in Eig.[^ where we plot the initial 
density profile along the x-axis for model ClBaObO of our 
sample (see Table |I] below). The AH is initially located at 
X = 0.46 and there exists a smooth density profile across 
the horizon. The remaining hydrodynamical quantities 
are treated in the same way. By using this approach we 
are able to evolve matter fields in the presence of punc¬ 
tures for very long times without the need for moving 
hydrodynamical excision zones. 


B. Tilted Kerr spacetime in improved 
quasi-isotropic coordinates 

We set up a Kerr black hole tilted about the x-axis by 
an angle /3o in the improved quasi-isotropic coordinates 
proposed by [68]. In those coordinates, the radius of the 
horizon does not shrink to zero in the extreme Kerr limit, 
but approaches r = M/4 > 0. The initial mass M of 
the Kerr BH is chosen to be equal to the value used in 
the ID calculation of the self-gravitating torus. The spin 
parameter a varies for different runs. The list of 21 models 
we use in our investigation is summarized in Table [IJ 
We perform the rotation of the BH about the x-axis 
in the following way: we first rotate by an angle /3o the 
coordinates x, 2 : to the tilted coordinates, 

= (^0 cos/3o -sin/?o) fy) . (21) 

yi j yO sin cos jd^ j J 
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TABLE I: Main characteristics of our initial models. From left to right the columns indicate the name of the model, the BH 
mass Mbh, the disk-to-BH mass ratio defined as the ratio of the irreducible mass of the central Schwarzschild BH and the 
total gravitational mass of the torus of the original initial data, the inner and outer disk radii rin and rout, the maximum 
rest-mass density pc, the polytropic constant of the EOS K, the orbital period Porb in code units and (milliseconds) and orbital 
frequency /orb at the radius of the initial pc, the specihc angular momentum prohle in the equatorial plane of the disk I (in 
terms of the Schwarzschild radial coordinate P), the BH spin parameter a, and the initial tilt angle /3o. All models are evolved 
with a E-law ideal fluid EOS with E = 4/3. Models indicated with a * were simulated using half the canonical resolution of 
Ax = 0.02. See main text for details. 


Name 

Mbh 

Q 

rin 

r out 

Pc 

(xio“b 

K 

Borb 

/orb 

(Hz) 

/-prohle 

a 

Pq 

(deg) 

D2a01b0 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.1 

0 

D2a01b5 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.1 

5 

D2a01bl5 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.1 

15 

D2a01b30 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.1 

30 

D2a05b5 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.5 

5 

D2a05bl5 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.5 

15 

D2a05b30 

L0002 

0.044 

3.42 

30.1 

1.05 

0.323 

150 (0.74) 

1360 

3.75 (const) 

0.5 

30 

ClBaObO 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.0 

0 

ClBaOlbS 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.1 

5 

ClBaOlblB 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.1 

15 

ClBa01b30 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.1 

30 

ClBa03b5* 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.3 

5 

ClBa03bl5* 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.3 

15 

ClBa03b30* 

0.9569 

0.160 

3.03 

22.7 

5.91 

0.180 

157 (0.77) 

1300 

3.67 (const) 

0.3 

30 

NClaObO* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.0 

0 

NClaOlbS* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.1 

5 

NClaOlblS* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.1 

15 

NCla01b30* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.1 

30 

NCla03b5* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.3 

5 

NCla03bl5* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.3 

15 

NCla03b30* 

0.9775 

0.110 

3.60 

33.5 

1.69 

0.170 

242 (1.19) 

843 

3.04 

0.3 

30 


We then calculate the spatial metric and extrinsic 
curvature Kij by performing coordinate transformations 
of the respective expressions given in [68] (in spherical 
polar coordinates) to our tilted Cartesian coordinates xL 
We set the initial shift to 0 and choose the following 
initial lapse profile: 


We find that this initial lapse profile greatly helps with 
the conservation of irreducible mass and spin during the 
initial gauge transition the system undergoes. The bigger 
the difference between the initially chosen lapse profile 
and the profile the system attains after the transition to 
the puncture coordinates, the larger the errors in irre¬ 
ducible mass and spin. As all quantities have been set up 
in the tilted coordinates x\ we also need to apply rota¬ 
tions to vectors and tensors, in order to obtain the correct 
components of the shift, spatial metric and extrinsic cur¬ 
vature in the original coordinate system that is going 
to be used in the simulations. We obtain the components 
of the shift vector with 

/3 = R^, (23) 

and the components of the spatial metric and extrinsic 


curvature with 

7 = Ro'R'^, (24) 

K = RKRT, (25) 

where R is the rotation matrix rotating from x* to xL 
Our choice of the initial BH spin magnitude is currently 
limited by our computational method and is therefore 
somewhat smaller than spins expected in certain astro- 
physical scenarios. For example, as the simulations of 
misaligned BH-NS mergers of [25H27] have shown, ac¬ 
cretion disks form in these systems only for relatively 
high (> 0.7) initial BH spins; see also [28| for disk mass 
predictions from BH-NS mergers. These initial spins in 
astrophysical merger scenarios are larger than the spins 
investigated in our study a G [0,0.5]. The reasons for 
our choice of smaller spins are the very demanding reso¬ 
lution requirements when evolving highly spinning BHs, 
see, e.g. [69]. Furthermore, spin conservation seems to 
be affected non-linearly with higher spin magnitude CO], 
which we have also observed in CD- The choice of astro- 
physically more realistic BH spin magnitudes would have 
been prohibitively expensive (see, e.g. [69] for details) and 
is beyond the scope of our study of self-gravitating ac¬ 
cretion disks around tilted Kerr BHs. Nevertheless, our 
results with dimensionless spin values of 0.5 should pro- 









FIG. 3: Contour plot in the xz plane of the initial rest-mass 
density profile in the disk for model ClBaObO. The innermost 
contour (blue line) corresponds to a density of 2.1 x 
while the outermost contour (red line) corresponds to 5.9 x 
10“^° (code units). The AMR grid structure, showing the 6 
innermost mesh refinement levels (out of 13), is shown as well 
in the figure. The box has a size of 46. See main text for details 
on the actual resolution of the different grids. 


vide a first qualitative understanding of what may hap¬ 
pen at dimensionless spin values of 0.7. We also note that 
the existence of the ’mass gap’ in BH masses in compact 
mergers has been questioned in m- While the peak of 
BH masses in these new considerations is still around 
SMq, very small BH masses in mergers are at least not 
ruled out as in previous studies. Furthermore, the BH 
mass predictions population synthesis models are based 
on assumptions about the exact details of the supernova 
explosion mechanism, see for instance m , where no mass 
gap is found using a delayed supernova explosion model 
assuming a 100-150 ms instability growth time. As we 
explained in the introduction, the smaller the BH mass, 
the smaller the BH spin needs to be in order to tidally 
disrupt the NS during the merger in order to leave a thick 
accretion torus around the central BH. It therefore seems 
that some of the combination of BH spins and mass-ratios 
studied in this work (while having been primarily moti¬ 
vated on computational grounds) are a possible (although 
maybe not the most likely) outcome from BH-NS merg¬ 
ers. Based on these considerations, model D2, the lightest 
in our study, is the most astrophysically relevant model of 
the three models studied. We choose the initial tilt angle 
/3o G [5°, 30°]. In the simulations of tilted BH-NS merg¬ 
ers of m, the authors found that for initial inclinations 
< 40°, the mass of the resulting accretion disk would be 
around 10 — 15% of the initial neutron star mass, whereas 
higher initial inclinations produced a sharp drop in the 
mass of the resulting disk. For those those pre-merger 


inclinations < 40°, the resulting tilt of the post-merger 
disk was around 10°. In the recent simulations of tilted 
BH-NS mergers of [27], the authors report a post-merger 
disk with a tilt ^ 20 — 30°, for an initial tilt of « 60°. 
These findings seem to be in accordance with the proba¬ 
bility distributions of post-merger tilt angles for BH-NS 
mergers obtained in m, which show a sharp cut-off in 
the probability distribution for post-merger tilt angles 
> 50°. 


C. Numerical setup 

The use of mesh-refinement techniques is of fundamen¬ 
tal importance in our simulations where different physical 
scales need to be accurately resolved. For this reason, we 
use the Carpet driver that implements a vertex-centered 
AMR scheme adopting nested grids m- 

Figure shows isocontours of the initial rest-mass den¬ 
sity profile for a representative disk model (ClBaObO) in a 
vertical cut (xz-plane) together with the initial innermost 
AMR grid structure. In our simulations, we use 13 levels 
of mesh refinement, with the outer boundary placed at 
4096. We perform a large series of simulations using two 
different resolutions, high resolution runs with a high¬ 
est resolution of Ax = 0.02 and simulations with half 
this resolution, see Table |l] Every refinement level ex¬ 
cept the innermost one is half the size of the preceding 
level. The innermost refinement level, as can be seen in 
Fig-i is more than half the size of the next refinement 
level. All simulations have been performed using a CFL 
factor of 0.25. In the outermost 6 refinement levels, we 
use the same time step as in level 7, in order to prevent 
instabilities in the spatially varying damping parameter 
T]. Using the same (smaller) time step in the outer levels 
is necessary, even though the spatially varying r] param¬ 
eter of m takes into account the time step limitation 
of the T] parameter described in ITT]. Details of the ac¬ 
curacy and convergence properties of our simulations are 
provided in Appendix 0 

IV. RESULTS 
A. Surface plots 

We start describing the morphology and dynamics of 
a representative model of our sample. The evolution of 
model ClBaOlbSO is shown in Fig.[^ This figure displays 
surface plots of the logarithm of the rest-mass density 
at six different stages of the evolution. The final panel 
shows the morphology of the system after the disk has 
completed 20 orbits around the central BH. The interpo¬ 
lation of the ID from QI spherical polar coordinates to 
Cartesian coordinates, along with the inclusion of a mod¬ 
erate BH spin in the originally non-rotating BH, results 
in a significant perturbation of the equilibrium model 
that triggers a phase of oscillations of the torus around 
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FIG. 4: Surface plots of the logarithm of the rest-mass density 
for model ClBaOlbSO at six different snapshots of the evolu¬ 
tion, t/torh = 0,4,8,12,16 and 20. The domain shown in all 
panels is 60 across the y-axis and 45 across the z-axis. The 
surface of the disk is set at a density of 10“^pc and half the 
disk is removed for visualization purposes. The BH lies at the 
center of each panel. The evolution shows the development 
of two accretion streams on to the BH as the disk inclina¬ 
tion increases during the evolution with respect to the initial 
equatorial plane. 


its equilibrium. These oscillations are present through¬ 
out the simulation and, as the disk is initially filling its 
Roche lobe, they induce a small accretion process of mat¬ 
ter through the cusp towards the BH. This does not re¬ 
duce significantly the total rest-mass of the torus as we 
show below. 

The oscillations of the disk are damped due to numer¬ 
ical viscosity and also by the formation of shock waves 
that propagate inside the disk and convert kinetic energy 
into thermal energy. These outward-propagating shocks, 
also found in the simulations of [18], are produced by 
in-falling material along the funnel walls which reaches 
inner high-density regions and bounces back. One such 
shock wave can be clearly seen on the right portion of 
the disk in the 12 orbits panel of Fig.|^ The initial tilt of 
the BH spin = 30°) twists and warps the innermost 
parts of the disk closer to the BH. This gradual effect 
becomes evident on the time series shown in Fig. [^and 
by the final 20 orbits the disk is significantly misaligned 
with the (horizontal) ^-axis. 

In Fig. we show surface plots of the rest-mass den¬ 
sity for the models GIB, NCI and D2 for the lower spin 


simulations (a = 0.1) at the final time of the evolution 
t/torh = 20. The effects of the initial BH tilt on the mor¬ 
phology of the disk at the final time are more clearly 
pronounced for the heavier models CIB (top row) and 
NCI (middle row). Although the spin of the BH is very 
low and the spacetime is therefore nearly spherically sym¬ 
metric, the initial tilt angle nevertheless causes a different 
evolution. In contrast, the significantly lighter model D2 
is seen to be hardly affected by the initial tilt by the end 
of the evolution, at least for these simulations with low 
initial spin. Similarly, in Fig. we plot the same rest- 
mass density surface plots as in Fig. for all our models 
with higher initial spin. Comparing the evolution of the 
disk morphology with the lower initial spin figure, we see 
that there are now significant changes in the morphology 
when increasing the initial tilt angle /3o for all models. 
The effects of evolving the disk in the tilted Kerr space- 
time are now more pronounced, as the higher spin causes 
a significant deviation from spherical symmetry in the 
spacetime. The growth of non-axisymmetric modes asso¬ 
ciated with the PPI (see below) in the two more massive 
models GIB and NCI causes an alignment of the overall 
disk angular momentum vector with the BH spin. This 
interesting effect will be discussed in Section [TV F| below, 
where we comment on the BH precession and nutation 
properties. 


In Figures and we plot linearly spaced isocon¬ 
tours and the logarithm of the rest-mass density in the 
plane perpendicular to the total angular momentum vec¬ 
tor of the disk. This kind of “equatorial” plots allows 
for a better visualization of the possible growth of non- 
axisymmetric structures in the disk. We note that due 
to the twisting and warping of the disk as a result of 
the initial BH tilt, the choice of the plane on to which 
project the isocontours for an adequate visualization is 
not straightforward and careless choices (as e.g. the sim¬ 
ple choice of the equatorial x^-plane) may hide impor¬ 
tant pieces of information on the dynamics and mor¬ 
phology. The top row of Fig. corresponds to mod¬ 


els ClBaOlbB (left), ClBaOlblB (middle), and ClBaOlbSO 
(right), and all of them are displayed at the time at which 
the growth of the fastest growing nona xisym metric struc¬ 
ture in the disk saturates (see Section IV D below where 
we discuss the mode growth of the PPI for the different 
models of our sample). The middle row of Fig. [^corre¬ 
sponds to models NClaOlbB (left), NClaOlblB (middle), 
and NClaOlbSO (right), also shown at the time the fastest 
growing mode has maximum amplitude, while the bot¬ 
tom row shows models D2a01b5 (left), D2a01bl5 (mid¬ 
dle), and D2a01b30 (right), which are displayed at the 
final time of the evolution (t = 20torb) when the ampli¬ 
tudes of the corresponding PPI modes (if present) are 
largest. Correspondingly Fig. shows the same type of 
isocontour plots as Fig. but for the models with higher 
initial spin. The top row of Fig. [^corresponds to mod¬ 


els ClBaOSbB (left), ClBaOSblB (middle), and ClBaOSbSO 
(right), the middle row of Fig. [^ corresponds to mod¬ 
els NClaOSbB (left), NClaOSblB (middle), and NClaOSbSO 
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FIG. 5: Surface plots of the (normalized) rest-mass density at the hnal time of the evolution t/torb = 20 for models GIB with 
a = 0.1 (top row), NCI with a = 0.1 (middle row), and D2 with a = 0.1 (bottom row). From left to right the columns correspond 
to initial tilt angles /3o = 5°, 15°, and 30°, respectively. The domain shown in all panels is 20 across the y-axis and 15 across 
the z-axis. The color palette and corresponding normalized density used are the same as in Fig. ^ 


(right) and the bottom row shows models D2a05b5 (left), 
D2a05bl5 (middle), and D2a05b30 (right). As in Fig. 
the times shown are at either saturation of the fastest 
growing non-axisymmetric structure or when the maxi¬ 
mum mode amplitudes are largest. 

The inspection of the different panels displayed in 
Figs. andshows noticeable differences in the flow mor¬ 
phology. On the one hand all six models GIB (top rows in 
both figures), corresponding to initial BH spins of a = 0.1 
and a = 0.3, respectively, exhibit a very prominent spi¬ 
ral arm once the stationary accretion phase is reached. 
This structure, which is visible for all three tilt angles 
considered, P = 5°, 15° and 30°, is associated with the 
growth of the m = 1 non-axisymmetric PPI mode, which 
is the fastest-growing mode and has the largest ampli¬ 
tude (see also Fig. 11). The middle rows, corresponding 
to the six models NCI, show the development of a spi¬ 
ral arm as well, however not as clearly pronounced as in 
the case of CIB. On the other hand the isodensity con¬ 
tours of all three models D2 displayed in the bottom row 
of Fig. for a BH spin a = 0.1 barely show any mor¬ 
phological deviation from the initial axisymmetry. Only 
the innermost regions closest to the BH show some slight 
non-axisymmetric structure, particularly for the model 
with the largest tilt angle {p = 30°) shown on the right 
panel. Model D2 is considerable lighter than model CIB, 


and model NCI has an initial non-constant specific angu¬ 
lar momentum profile and according to m more massive 
tori with j = const profiles favor the appearance of the 
PPI with respect to less massive tori (see in particular 
Figs. 2 and 3(a) of [19] for their model Cl which has very 
similar properties regarding size and BH-to-disk mass ra¬ 
tio to our model CIB). Our results seem to only partially 
confirm those previous findings as there seems to be some 
weak dependence of a moderate m = 2 PPI growth with 
increasing tilt angles. 


The dependence on the BH spin seems however to be 
more significant for the light model D2, as shown in the 
panels at the bottom row of Figures and For the 
higher spin runs shown in Fig. all three panels show 
non-axisymmetric morphological features. All three snap¬ 
shots correspond to the same final time of the evolution 
{t = 20torb) and the dynamical differences between the 
three models are only due to the initial BH tilt angle. 
For ^ = 5° (left panel) the m = 1 spiral structure is 
visible and dominant but, at the same time, the m = 2 
mode seems to have reached an almost similar amplitude 
(see the bottom panel of Fig. 11). Such m = 2 feature 
becomes clearly dominant the larger the BH inclination 
becomes, as depicted in the central regions of the middle 
and right panels. 
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FIG. 6: Surface plots of the (normalized) rest-mass density at the hnal time of the evolution t/torb = 20 for models GIB with 
a = 0.3 (top row), NCI with a — 0.3 (middle row), and D2 with a — 0.5 (bottom row). From left to right the columns correspond 
to initial tilt angles /3o = 5°, 15°, and 30°, respectively. The domain shown in all panels is 20 across the y-axis and 15 across 
the z-axis. The color palette and corresponding normalized density used are the same as in Fig. ^ 


B. Maximum rest-mass density evolution 


We next show the evolution of the maximum rest-mass 
density in the disks in Fig. normalized to the corre¬ 
sponding initial central density of each model (see Ta¬ 
ble]^. For model D2 (bottom panel), the effect of tilting 
the BH has very little effect on the evolution of pmax for 
a BH spin of a = 0.1. However, for a = 0.5, the evolution 
of pmax differs for the different initial tilt angles. Model 
D2a05b5 exhibits an upward drift from 10 orbits onwards 
to bring pmax to the level of the lower spin runs towards 
the end of the simulation. This might be connected to 
the growth of the m=l non-axisymmetric mode in this 
model, as discussed in the previous section. In the initial 
stages of the evolution, pmax is higher for larger tilt an¬ 
gles, a dependence that is also observed in the accretion 
rate (see Section IVC). The maximum rest-mass density 
stays below its initial value for almost the entire evolu¬ 
tion for both models D2 and NCI (middle panel). This is 
very different in model CIB (top panel) where the occur¬ 
rence and saturation of the PPI has a drastic effect on 
the evolution of pmax- When the PPI starts growing, a 
situation which is accompanied by an exponential growth 
of the accretion rate, pmax grows to up to twice its initial 
value, reaching the highest value for the untilted model 


ClBaOO. It subsequently settles to a value that is very 
similar for all initial spins and tilt angles, but the time 
of growth and the shape of the “lump” in the density 
evolution are different for the different spins and initial 
tilt angles. Model NCI shows features from both mod¬ 
els D2 and CIB described above. Similarly to model D2, 
Pmax stays below its initial value for almost the entire 
evolution, with model NClaOO being an exception where 
it grows above the initial value briefly. There is, as in 
models CIB, a rapid growth in the central density, albeit 
somewhat smaller for the non-constant angular momen¬ 
tum disks, associated with the growth of the PPI in all 
the models. The magnitudes at which the mode growth 
saturates are similar, while the time at which the growth 
sets in depends on the initial spin and tilt angle. 


C. Evolution of the accretion rate 

Fig, pr] shows the time evolution of the mass accretion 
rate for all models of our sample, in units of Mq/s. The 
colors used for the different models as well as the line 
style follow the same convention as in the rest-mass den¬ 
sity evolution figure (Fig. |^. The top panel displays the 
mass flux for models CIB, the middle panel for models 
NCI, and the bottom panel for models D2. The mass flux 
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FIG. 7: Linearly spaced isocontours in plots showing the logarithm of the rest-mass density in the plane perpendicular to the 
total angular momentum vector of the disk. From left to right the rows show models ClBaOlbS, ClBaOlblS and ClBaOlbSO 
(top), NClaOlbB, NClaOlblb and NClaOlbSO (middle), and D2a01b5, D2a01bl5, and D2a01b30 (bottom). The color palette and 
corresponding normalized density used are the same as in Fig. The morphological structures shown in the different panels 
reflect the dominant PPI mode at the time of maximum non-axisymmetric mode amplitude. See main text for further details. 


is computed as the instantaneous flux of matter across 
the AH surface which is parameterized in spherical coor¬ 
dinates by the polar angle 0 and azimuthal angle 0 and 
is given by: 

^277 nTT 

M = 27rr^ / / D d(j)dO ^ (26) 

io Jo 

where r is the radius of the AH and is the radial 
velocity of the fluid crossing the sphere. All three panels 
of Fig. show that after a transient initial phase the 
mass accretion rate is seen to tend asymptotically to a 
fairly constant value. 

In the previous fixed spacetime simulations of [32] it 
was found that the tilt and the twist of the disks are 
not strongly responsive to M. Our results for a dynami¬ 
cal spacetime validate those earlier findings, as the final 
dispersion found in M once this quantity reaches steady- 


state values is relatively small irrespective of the BH spin 
and inclination angles (M [Mq/s] G [0.2,1.5] for models 
CIB and NCI, and [0.06,0.2] for model D2). However, the 
transition to the steady-state does seem to show a clear 
dependence on the BH spin. We find that the initial M 
values are consistently smaller the larger the BH spin, 
about 2 orders of magnitude for models D2 after or¬ 
bit and about half that value for the more massive disks 
CIB and NCI. Another trend, which is most clearly seen 
in model D2, is the initial dependence of the accretion 
rate on the tilt angle as well. We find that the higher the 
initial tilt angle, the higher the accretion rate in the early 
stages of the evolution. This is because the innermost sta¬ 
ble circular orbit (ISCO) of a rotating Kerr BH has an 
ellipsoidal shape, attaining the largest size at the equato¬ 
rial plane of the BH. For prograde orbits such as those in 
our work, the higher the BHs spin, the smaller the size of 
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FIG. 8: Linearly spaced isocontours in plots showing the logarithm of the rest-mass density in the plane perpendicular to the 
total angular momentum vector of the disk. From left to right the rows show models ClBa03b5, ClBa03bl5 and ClBa03b30 
(top), NCla03b5, NCla03bl5 and NCla03b30 (middle), and D2a05b5, D2a05bl5, and D2a05b30 (bottom). The color palette and 
corresponding normalized density used are the same as in Fig. The morphological structures shown in the different panels 
reflect the dominant PPI mode at the time of maximum non-axisymmetric mode amplitude. See main text for further details. 


the ISCO and the closer can the matter orbit around the 
BH before falling in. When tilted, this centrifugal barrier 
becomes however less and less effective and the accretion 
rate becomes therefore higher for larger inclinations. 


A common feature present in the mass-flux plots of 
the simulations of [32] was a late-time “bump” which the 
authors attributed to the reflection of an outgoing wave 
within the disk at the outer edge. That feature does not 
seem to be so apparent in our models apart from maybe 
model CIB which shows a bump around t ^ lOtorb some¬ 
what similar to that reported by [32]. However, in our 


interpretation, the phase of rapid growth in M in less 
than about 5 orbits seems to have a different origin. This 
conclusion is reached by comparing the time at which the 


growth of the PPI saturates (discussed in Figs. 11 and 12 


in the following section) with the saturation time of the 
accretion rate for model CIB, which yields a significant 


close agreement. We note that this peak in the accretion 
rate is the highest of all our models, ^ BOMq/s, lasting 
for less than about 1 ms. Therefore, we identify the ex¬ 
ponential growth and subsequent saturation of the mass 
accretion rate with the growth and saturation of the PPI 
in our models. 


The total rest mass accreted during the evolution 
ranges within Mace ^ [0.022,0.033], [0.005,0.013] and 
[0.0006,0.002] for models CIB, NCI and D2, respectively. 
The accretion is facilitated by the outward transport of 


angular momentum, which we describe in section IV H 


below. The growth of the PPI, to which we attribute the 
high peak accretion rates, is very effective in transport¬ 
ing angular momentum outwards. While we have per¬ 
formed our study without considering magnetic fields, it 
is known that magnetic fields are very effective at trans¬ 
porting angular momentum in accretion disks m , via 
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FIG. 9: Evolution of the maximum rest-mass density, normal¬ 
ized by its initial value, for models CIB (top), NCI (middle) 
and D2 (bottom). The dashed curves correspond to the high 
spin models of our sample. 


the development of the so-called magnetorotational in¬ 
stability [73110]. Not taking into account the effects of 
magnetic fields in accretion disks seems to underestimate 
the accretion rate [81], see also recent GRMHD simula¬ 
tions of BH-NS mergers (82] [83] and the accretion rates 
reported therein. 


D. PPI growth in the disk and BH movement 

In their seminal paper on the PPI m, the authors 
showed that vertically thick, radially slender tori with 
constant specific angular momentum profiles are unsta¬ 
ble against the development of global non-axisymmetric 


FIG. 10: Evolution of the accretion rate in Mq/s for models 
CIB (top), NCI (middle) and D2 (bottom). 


modes. Two of our models, CIB and D2, have constant 
specific angular momentum profiles, while model NCI 
has a radial dependence of / ^ where q = 0.11. 
In [84], the authors showed that tori with non-constant 
angular momentum profiles might become PP-unstable 
as well. In order to check for the onset and the growth 
of non-axisymmetric instabilities, we monitor the mode 
amplitudes computed using Eq. (36), which is basically 
a Fourier decomposition of the azimuthal distribution of 
the density in the disk [84]. In simulations of accretion 
tori susceptible to the development of the PPI performed 
on spherical grids, it is customary to induce the instabil¬ 
ity with small initial non-axisymmetric density or veloc¬ 
ity perturbations. As outlined any non-axisymmetric per¬ 
turbation will trigger the growth of the instability pro¬ 
vided the torus is not stable against it. We do not ac- 
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FIG. 11: Evolution of the first four non-axisymmetric (PPI) 
modes scaled by the total rest-mass density (Dq) for mod¬ 
els ClBaOlbSO (top), NClaOlbSO (middle) and D2a01b30 (bot¬ 
tom). The inset shows the evolution of the m = 1 mode and 
the polar distance of the centre of the BH from the origin, 
both having been rescaled with their respective maximum val¬ 
ues attained during the evolution. 


tively seed the PPI in any of our models. Instead, there 
are two main sources of perturbations in our simulations, 
both connected to the fact that we evolve the BH-disk 
system in a Cartesian grid. First, the exact axisymme- 
try of the initial data is lost when we interpolate the 
initial data onto the Cartesian grid of the evolution, in¬ 
troducing small, non-axisymmetric perturbations in all 
interpolated quantities. This is due to the fact that in 
Cartesian grids neighboring cells lie on concentric shells 
only in the limit of infinite resolution. The second source 
of perturbation is the so-called junk-radiation leaving 
the BH while the gauge evolves from the values speci¬ 


fied by the initial data to the puncture gauge attained 
during the evolution. Here, we need to distinguish be¬ 
tween models evolved around Schwarzschild BHs and the 
models evolved around tilted Kerr BHs. For the former, 
the junk-radiation leaving the BH initially should be ax- 
isymmetric. Indeed, when analyzing this initial burst of 
radiation with the Weyl scalar ^ 4 , we find that the am¬ 
plitudes of the ^ 4 ^ multipoles (which are axisymmetric) 
are the largest. Nevertheless, the initial junk-radiation 
has also non-vanishing amplitudes in the other ^ 4 "^ mul¬ 
tipoles. Again, the reason for the presence of these non- 
axisymmetric ^ 4 ^^ multipoles seems to be connected to 
the evolution in the Cartesian grid, and specifically to the 
existence of mesh refinement boundaries, which are not 
axisymmetric even in the continuum limit (see [85] for 
interference patterns produced by the radiation crossing 
mesh refinement boundaries). Finally, in the case of the 
disk being around a tilted Kerr BH, the junk-radiation 
reaching the disk is now truly non-axisymmetric even 
without taking into account the effects of the Cartesian 
grid and the mesh refinement boundaries, as the axis of 
rotation is tilted with respect to the z = 0 plane. 

We show in Fig. [^the evolution of the amplitude of 
the first four non-axisymmetric modes Dm (m = 1 — 4) in 
the disk for models ClBaOlbSO (top), NClaOlbSO (mid¬ 
dle) and D2a01b30 (bottom). Those amplitudes are com¬ 
puted using Eq. Q. In addition, in the inset of all three 
plots we also show as a black solid curve the evolution of 
the polar distance of the center of the BH from the origin 
of the computational grid, which we label ppH- 

Model ClBaOlbSO is the one that most clearly is found 
to develop the PPI, as signaled by the exponential growth 
of the m = 1 mode. Modes m = 2 — 4 also show a phase 
of exponential growth lasting for about 2 orbital periods 
(between t ^ 8 — lOtorb) but the saturation amplitudes 
reached are 1 to 2 orders of magnitude smaller than for 
m = 1 . As mentioned in the previous section, the growth 
of the PPI is accompanied by an exponential growth of 
the mass accretion rate (see Fig. 10). For this model the 
PPI clearly saturates at around t ^ lOtorb? when the 
amplitude of the m = 1 mode quickly drops by an order 
of magnitude. After that, the amplitude of all the modes 
we have analyzed remains rather constant. The strength 
of the m = 2 — 4 modes becomes similar at late times 
while the m = 1 mode still retains a significantly larger 
amplitude. 

The polar distance of the BH from the origin for model 
ClBaOlbSO closely follows the behavior of the m = 1 
fastest growing PPI mode. This distance is found to grow 
at the same actual rate as the m = 1 mode, because the 
circular motion of the m = 1 overdensity hump (well 
visible in the top-rightmost panel of Fig. causes the 
BH to start moving in a spiral trajectory. After mode 
saturation the polar distance of the BH from the origin 
remains almost constant, reflecting a similar trend in the 
m = I evolution. This spiral motion of the BH as a result 
of the development of the PPI has also been observed in 
the simulations of m- 
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FIG. 12: Evolution of the non-axisymmetric m = 1 mode for 
all models, namely CIB (top), NCI (middle) and D2 (bottom), 
rescaled to the maximum amplitude of the Di mode of the 
respective untilted models. 


For the non-constant angular momentum model 
NClaOlbSO (middle panel) the growth rate of the m = 1 
mode is smaller, although it is still the dominant mode, 
as in model ClBaOlbSO. The growth of the mode satu¬ 
rates somewhat later at around t ^ 12.5 torb but, con¬ 
trary to the evolution of ClBaOlbSO, it does not drop 
significantly after saturation, remaining fairly constant 
thereafter. Furthermore, it remains the dominant mode 
at late times while the three remaining modes attain sim¬ 
ilar amplitudes. The evolution of the polar distance of the 
BH from the origin shows a small secular drift after the 
saturation of the PPL 

On the other hand, model D2a01b30 (bottom panel) 
shows a very distinct behavior. The amplitude of all four 


modes and the polar distance of the BH remain below 
the values attained in models ClBaOlbSO and NClaOlbSO, 
typically in the range between 1 to 2 orders of magnitude. 
Furthermore, the m = 1 and m = 2 modes grow at ap¬ 
proximately the same rate and to similar final values dur¬ 
ing the evolution. The growth rate of ppH is very similar 
to the growth rate of the modes m = 1 and m = 2. We 
also note that at early times the m = 4 mode shows the 
largest amplitude (purple line; not so clearly noticeable 
in the previous two models), perhaps an artifact of the 
Cartesian grid we use in our simulations. Nevertheless, 
the amplitude of this mode at late times is sufficiently 
smaller than the dominant amplitudes (mostly that of 
the m = 1 mode) to safely consider its effect on the dy¬ 
namics negligible. From the mode analysis, we can there¬ 
fore say that model D2 is essentially PP-stable. We will 
return to this point in section [TVH[ where we analyze the 
angular momentum transport in the models. 

The results we have just described connect the expo¬ 
nential growth and saturation of the mass accretion rate 
with the growth and saturation of the PPI in the models 
where it has unambiguously taken place. To further jus¬ 
tify this, we plot in Fig. [^the evolution of the m = 1 
mode for all models considered in our sample. Focus¬ 
ing on the evolution of the models CIB displayed in this 
figure (top panel) we can clearly see that the temporal 
order of the saturation of the mode growth closely coin¬ 
cides with the various peaks present in the mass accretion 
rate evolutions shown in Fig. (compare curves of the 
same color). The comparison among the three types of 
models also shows that the saturation mode growth for 
models CIB is only slightly larger than for models NCI, 
their late-time values being fairly similar, but about 2 or¬ 
ders of magnitude larger than in models D2. The slope of 
the mode growth is also steeper in models CIB, while the 
behavior of the mode evolution in the case of the light 
tori D2 almost seems to indicate that the development of 
non-axisymmetric instabilities is only marginal for such 
models. A striking feature worth highlighting from the 
middle panel of Fig. is the non-existing dependence of 
the late time m = 1 mode amplitude on the BH spin and 
tilt angle for the non-constant angular momentum tori 
NCI. 


To demonstrate the connection between the growth of 
non-axisymmetric modes in the disk with the movement 
of the central BH, we plot in Fig. [^the position of the 
BH center on the xy plane for models ClBaOlbSO (top), 
NClaOlbSO (middle) and D2a01bS0 (bottom). This fig¬ 
ure shows several interesting features connected to the 
growth of the modes in the tori. The motion of the BH 
is caused by the formation of a “quasi-binary system” 
between the BH and the m = 1 overdensity lump. It is 
this binary motion that causes the long-term emission 
of GWs reported in m- We will return to this in the 
discussion of the GW emission connected to the PPI in 
section |IVT| For model ClBaOlbSO, the development of 
the dominant m = 1 structure in the disk exerts a small 
kick in the BH. As a result the BH starts moving in a 
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FIG. 13: Position of the BH center on the xy plane for models ClBaOlbSO (left), NClaOlbSO (middle) and D2a01b30 (right). 
For models ClBaOlbSO and NClaOlbSO, the red dot indicates the point in the trajectory when the PPI saturates. Note that the 
scale ratio is not preserved in the hgures for clarity in the visualization and that the length scale is signihcantly different in 
every panel. 


spiral trajectory until the PPI saturates and the mode 
amplitude drops, at which point the BH is roughly lo¬ 
cated at {x^y) = (—0.3,—0.5). After saturation, the xy 
position of the BH is a combination of linear and cir¬ 
cular motion. In model NClaOlbSO, the linear motion is 
stronger from the beginning, causing a linear shift in the 
growing spiral motion. Upon saturation of the PPI, the 
circular motion continues with a smaller radius, in accor¬ 
dance to the evolution of the m = 1 mode which remains 
at a rather constant value after saturation and up until 
the end of the evolution. Finally, for model D2a01b30 we 
have a superposition of circular and linear motion again, 
as the modes m = 1 and m = 2 are of almost equal mag¬ 
nitude during the evolution. Here, the overall distance 
covered by the BH is significantly smaller than for the 
heavier models ClBaOlbSO and NClaOlbSO, as the mode 
amplitudes are much smaller for D2 and the disk is also 
lighter. 

We argue that the non-axisymmetric instability we ob¬ 
serve in our models is indeed the PPI based on the az¬ 
imuthal mode evolution in the disk and the development 
of the m = 1 over density lump. However, we do not claim 
to have presented a direct, undeniable proof that the in¬ 
stability observed in our simulations is the analytical PPI 
studied in m- We also note that there exist other insta¬ 
bilities in self-gravitating disks which display spiral arm 
formation in the disk or even disk fragmentation, such as 
the axisymmetric local Toomre instability [86] , but as m 
have shown, the Toomre instability is unlikely to affect 
radially slender accretion tori such as the ones considered 
in our study. 


E. PPI saturation and black hole kick 

The complete 3D trajectory of the BH and its spin vec¬ 
tor is plotted for model ClBaOlbSO in Fig.[^ In this plot 


each small sphere with an attached arrow corresponds to 
the position of the BH and the magnitude (scaled for bet¬ 
ter visualization) and direction of its spin vector. The po¬ 
sitions are plotted at equal time intervals, so we see that 
the BH is moving much faster during the final stages of 
the PPI than before its occurrence or after its saturation. 
We can also see that there is a movement in the vertical 
direction, which is caused by the fact that the “binary” 
motion of the BH and m = 1 overdensity lump does not 
lie in the equatorial plane anymore due to the initial tilt. 
While the PPI is growing, the BH-torus system therefore 
moves up and down in an oscillatory fashion. As we have 
seen in the m = 1 mode plots for models CIB (see the top 
panel of fig. 12), the saturation of the PPI is very fast. 
The rapid saturation of the PPI (and the corresponding 
destruction of the m = 1 overdensity lump) result in a 
mild vertical kick of the BH-l-torus system. To see that 
this is indeed connected to the saturation of the PPI, we 
plot in Fig. the linear momentum radiated away by 
gravitational waves in the ^-direction for models ClBaOl. 
From the plot, we can clearly see that the time of max¬ 
imum emission of linear momentum corresponds exactly 
to the time the PPI saturates for each model. The ra¬ 
diated linear momentum has been calculated with the 
pyGWAnalysis package [88] . 


F. BH precession and nutation 

Contrary to the test-fluid simulations of [32l [33], as 
we are evolving the full spacetime, we can monitor the 
response of the BH to the evolution of the disk in the ini¬ 
tially tilted configuration. In particular we can measure 
the total precession of the BH spin about the z-axis and 
its rate, as well as the inclination of the spin with respect 
to the 2 :-axis and the corresponding nutation. 

In Fig. I^we plot the evolution of the precession rate 
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FIG. 14: Full 3D trajectory of the BH along with its spin 
vectors for model ClBaOlbSO. The magnitude of the spin vec¬ 
tors has been rescaled with a constant factor for visualization 
purposes. 



FIG. 15: Radiated linear momentum in the z-direction for 
models ClBaOl. 


(Jbh of the central BH about the z-axis for models CIB 
(top), NCI (middle) and D2 (bottom panel) in radians 
per orbital period. The first thing to mention is that the 
constant angular momentum disks D2 and CIB show a 
clear dependence of the precession rate with the initial 
BH spin magnitude. Namely, for models D2, the larger 
the spin the smaller the precession rate, and for models 
CIB smaller spins lead to higher precession rates. How¬ 
ever, for the non-constant angular momentum torus NCI 
the precession rate is very similar irrespective of the spin 
magnitude, particularly at the early stages (t < lOtorb)- 
In models D2 the evolution of (Jbh is very smooth, and 
changes sign for the higher spin runs after about 17 or¬ 
bits, whereas it remains fairly constant for the lower spin 
simulations. 

A fact worth emphasizing is that there is essentially 
no dependence of the precession rate on the tilt angle 
for constant angular momentum tori. This is particularly 
true for models D2 but also models CIB show this lack of 
dependence up until the growth of the PPL In all models 



FIG. 16: Evolution of the precession rate of the BH spin about 
the z-axis for models CIB (top), NCI (middle) and D2 (bottom) 
in units of radians/Grb • 


CIB we observe a prominent modulation of the precession 
rate when the PPI enters the final stages of its growth 
and then saturates, as well as the change in sign for the 
higher spin runs. The magnitude of the precession rate is 
about an order of magnitude higher compared to mod¬ 
els D2, which we attribute to the higher disk-to-BH mass 
ratio models CIB possess. The evolution of the BH pre¬ 
cession rate for models NCI seems to fit in between the 
two other models, the magnitude of the precession rate is 
somewhere in between but not very different evolutions 
are observed for the two BH spins considered. We also 
note from Fig. [^that models ClBaOl and D2a01 attain 
a fairly constant precession rate by the end of our sim¬ 
ulations, namely ^ 20Hz and ^ 6Hz for models ClBaOl 
and D2a01, respectively. 
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FIG. 17: Total precession of the BH spin about the z-axis for 
models CIB (top), NCI (middle) and D2 (bottom) in units of 
radians. 
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FIG. 18: Evolution of the nutation rate of the BH spin about 
the z-axis for models CIB (top), NCI (middle) and D2 (bottom) 
in units of radians/Grb • 


The total precession of the BH spin about the z-axis 
is shown in Fig. The lower spin models ClBaOl, that 
show the highest BH precession rate, have completed a 
quarter of a full precession cycle by the end of our sim¬ 
ulations at t = 20torb- Given that the precession rate for 
those models appears to settle to a rather constant value 
(see Fig. 16) we would need a prohibitively expensive 
evolution of about 80 orbital periods for a full precession 
cycle of the central BH for those models. If we assume 
that the precessing Kerr BH radiates GW as a freely 
precessing rigid body, it will radiate at dpH and 2 (Jbh 
in the / = 2, m = 1, 2 multipole modes [89]. This means 
we would have to evolve the BH-torus system for at least 
one full BH precession cycle in order to see the slow mod¬ 
ulation of the GW signal caused by the BH precession. 
It remains an interesting open issue to see if the upcom¬ 


ing advanced GW detectors will be able to measure these 
GWs in the low tens of Hz. 

The imprint of the growth of non-axisymmetric insta¬ 
bilities in the disk on the BH response can be seen even 
clearer in Fig. This figure shows the evolution of the 
nutation rate for all models, that is, the evolution of the 
temporal variation of the tilt angle of the BH with the 
z-axis, z>bh- Foi* models CIB (top panel) there is a clear 
period of alignment of the BH spin with the z-axis around 
the time when the PPI starts to grow. This is signaled 
by the rapid fall and rise of the nutation rate around 
the lOtorb mark for the slow spin (a = 0.1) BHs, but 
it is also (less) visible in the more rapidly rotating BHs 
(dashed lines). Upon saturation of the instability the nu¬ 
tation rate becomes almost constant and fairly close to 
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FIG. 19: Evolution of the tilt angle of the BH spin with re¬ 
spect to the z-axis for models CIB (top), NCI (middle) and D2 
(bottom) in units of radians. 


being again zero for the a = 0.1 BH models (see solid 
lines). This behavior is consistent with the evolution of 
the BH tilt angle itself displayed in Fig.[^ which, for the 
case of models CIB with small spins, shows two distinct 
zones of constant spin inclination, particularly for models 
with Pq = 5° and 15° (model with /3o = 30° still shows 
a negative slope at late times). The transition to zero 
nutation rate in the a = 0.3 CIB models after PPI satu¬ 
ration seems to take a longer time than that displayed in 
Fig.[T8l 

In addition Fig. [T^ also shows that the maximum nu- 
tation rate attained increases with increasing initial tilt 
angle. This is very clearly seen in models CIB but it is also 
visible in most of the models D2 and NCI. In particular the 
evolution of z>bh for the non-constant angular momentum 


models NCI with a = 0.3 (dashed curves) does not seem 
to reach a local maximum, at least in the timescales we 
can afford in our simulations, despite these models are 
also affected by the PPI. This might be explained by the 
fact that the non-axisymmetric m = 1 mode does not 
drop as sharply after the saturation of the PPI in the 
models NCI as it does in the CIB models (see Fig. 12). 
The growth of the mode in the disk seems to cause an 
alignment of the BH spin with the 2:-axis. Furthermore, 
the nutation rate is much larger for the CIB models. Fi¬ 
nally, models D2, in the absence of any strong growth of 
PPI modes, do not show strong modulation of the nuta¬ 
tion rate. The modulation of this rate is therefore closely 
connected to the existence of m = 1 non-axisymmetric 
modes in the disk. 

The corresponding evolution of the inclination of the 
BH spin from the z-axis is seen in Fig. The panels 
show that a significant realignment of the BH spin and 
the 2:-axis takes place by the end of the simulations in 
models CIB, i.e. in those that develop the PPI most sig¬ 
nificantly. For all models we see that the final inclination 
is smaller the bigger the initial tilt angle is, for equal 
spin magnitudes. This trend also seems to be more pro¬ 
nounced the higher the initial spin magnitude is. 


G. Disk twist and tilt 


We turn now to describe in this section the response 
of the disks as they evolve in the tilted spacetime, using 


the diagnostics we have introduced in Section |HD[ In 
Figures and we plot the evolution of the tilt and 
precession of the total angular momentum vector of the 
disk Joisk- This vector is the sum of the angular momen¬ 
tum vectors of each shell, as described in Section |nD} 
We will refer to these as the global disk tilt and global 
disk precession around the BH spin, in order to distin¬ 
guish them from the analysis of the twist cr(r) and tilt 
z/(r) angles in the individual disk shells. 

The evolution of the global disk tilt, shown in Fig. 
shows resemblances with the evolution of the tilt of the 


BH spin (see Fig. 19). For models CIB (top panel) we see 
that the development of the PPI not only causes a partial 
realignment of the BH spin with the 2:-axis, but also a re¬ 
alignment of the BH spin and disk angular momentum. 
This alignment seems to become more constant towards 
the end of the simulation than the alignment of BH spin 
and the 2:-axis. We also note that the amount of realign¬ 
ment is reversed for the global disk tilt for models CIB. 
Here the higher initial spin magnitude leads to a lower 
amount of alignment, in contrast with the evolution of the 
BH tilt. On the other hand, the alignment between BH 
spin and total disk angular momentum for models NCI 
(middle panel) is not as pronounced as the realignment 
of the BH spin with the 2:-axis. As in the case of models 
CIB, the trend of more alignment with higher spin for the 
BH tilt is inverted for the evolution of the global disk tilt. 
There are also oscillations in the evolution of the global 
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FIG. 20: Evolution of the tilt angle z^Disk between then BH 
spin and the total disk angular momentum Joisk for models 
GIB (top), NCI (middle) and D2 (bottom) in units of radians. 


FIG. 21: Evolution of the total precession of the total disk 
angular momentum Joisk around the BH spin for models CIB 
(top), NCI (middle) and D2 (bottom) in units of radians. 


disk tilt for models NCI, which are not visible in the other 
two set of models. These oscillations seem to be caused 
by the persistent m = 1 non-axisymmetric mode in mod¬ 
els NCI. Finally, for models D2 (bottom panel), neither 
the initial spin magnitude nor the initial tilt angle seem 
to affect the evolution of the global disk tilt, which is 
furthermore hardly changed during the entire evolution 
for all initial spins. The development of the PPI there¬ 
fore seems to cause an alignment of the BH spin with the 
total disk angular momentum, which is larger for smaller 
initial spins. As both z^bh and z^Disk decrease during the 
growth of the PPI, it is clear that the BH spin is tilted 
into the orbital plane of the disk. 

In Fig.j^we plot the global precession of the total disk 
angular momentum vector about the BH spin axis for 


models CIB (top), NCI (middle) and D2 (bottom panel). 
For all models, the higher the spin (dashed curves), the 
larger the global precession, as expected. The initial tilt 
angle does not seem to influence strongly the evolution 
for models CIB and D2, while the evolution for NCI shows 
oscillations towards the end of the evolution. Note that 
these oscillations are superimposed on the slow growth 
of the precession, which always has a positive slope, as 
expected. The oscillations are caused by the wobbling 
and smaller precession cycles about the direction of the 
global vector. These oscillations are then visible in the 
projection of the disk angular momentum vector onto 
the equatorial plane of the BH. The amplitude of the 
oscillations is larger the smaller the initial tilt angle. The 
non-axisymmetric modes that survive much longer in the 
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FIG. 22: Spacetime t — r diagrams of the radial disk tilt profile i>{r) for models ClBaOl, ClBaOS (top two panels), NClaOl, NClaOS 
(middle two panels), and D2a01 and D2a05 (bottom two panels). From left to right the panels show the initial tilt angles /3o of 
5°, 15° and 30°. Each panel shows the radial profile of the disk tilt for the entire time evolution. 


case of models NCI could be causing these oscillations in 
the global disk tilt and precession. 


In Fig. 22 we show spacetime (t — r) diagrams of the 
radial disk tilt profile for the entire time evolution for all 


tilted models in order to see if and which disk models 
become warped during their evolution. The initial radial 
tilt profile is flat and the panels cover the initial radial ex¬ 
tent of each model. For models CIB (top two rows) we see 
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that by the end of the simulations the radial tilt profile for 
larger radii is significantly lowered compared to the initial 
tilt angle. Furthermore, there is a clear distinction in the 
radial tilt profile before and after the saturation of the 
PPI for the lower spin models (a = 0.1), in accordance 
with what we have observed for the global disk angular 
momentum tilt (see Fig. 20). For the higher spin models 
(a = 0.3) the tilt profile evolution is broken up into two 
regions as well. In these cases, the peak close to the ori¬ 
gin reemerges after the saturation of the PPI. Likewise, 
models NCI (two middle rows) also show a peak close to 
the origin, but contrary to the previous set of models, the 
late time tilt profiles are oscillating in the regions closest 
to the central BH. The oscillations are stronger for both, 
radii closer to the origin and smaller initial BH spin. As 
described before, we believe that these oscillations are a 
consequence of the long-lasting non-axisymmetric m = I 
mode in the disk. 


The oscillations of the tilt angle close to the origin in 
models NCI might be connected to the so called Kozai- 
Lidov (KL) mechanism [90l |9T], where the eccentricity 
and inclination of particle orbits around an object that 
is itself in a binary are interchanged in an oscillatory 
fashion. Recently [92] performed simulations of the KL 
mechanism in hydro dynamical disks around a component 
of a binary system, finding oscillations of the inclination 
angle. The persistent m = 1 structure in the disks of 
models NCI can be thought of as an eccentric distribution 
of matter. While our disks are not evolved around a BH 
in a binary system, we have nevertheless shown in Sec- 
tion |IV Pj that the growth of the m = 1 non-axisymmetric 
modes causes the BH to start moving in a “quasi-binary” 
with the overdensity lump in the disks. As seen in Fig.[T^ 
for models NCI, the long lasting “planet” causes the BH 
to move in a continuous elliptic fashion in the projection 
onto the x^-plane. We believe this movement of the BH 
induced by the PPI serves as the correct boundary con¬ 
ditions for the KL mechanism to occur in models NCI. 
In [92] the authors estimate the timescale of the KL os¬ 
cillations in their disks as being tkl ^ IbPb, where P 5 is 
the orbital period of the central object in the binary. In 
our case, the periods of the tilt angle oscillation and the 
orbital motion of the Kerr BH are roughly equal for all 
models. We plan to further investigate the correctness of 
this interpretation of the tilt angle oscillations in a future 
work. 

To complete the analysis of Fig.[^ we note that mod¬ 
els D2 (two bottom rows) become clearly warped, with 
all models showing a peak near the origin. The profiles 
of the low spin models D2a01 are very similar in shape 
for the three initial tilt angles. The peak in models D2a05 
is seen to oscillate stronger the smaller the initial tilt an¬ 
gle. Note that for none of the models we see an alignment 
of the inner region of the disk with the equatorial plane 
of the BH (this would mean z^(r) = 0 there). We there¬ 
fore find no occurrence of the Bardeen-Petterson effect 
in our simulations, at least in the timescales considered. 
This is in agreement with the results found in the fixed¬ 


spacetime inviscid simulations with no angular momen¬ 
tum transport of [32] and also with the GRMHD sim¬ 
ulations of [ 33 ] that included angular momentum trans¬ 
port driven by the magneto-rotational instability. The 
observed profile and evolution of the disk warp seems to 
be in qualitative agreement with the analytic work on 
warp propagation as bending-waves in tilted thick ac¬ 
cretion disks; see, for instance, the analytic tilt profiles 
in [93] , the analytic model of linear warp evolution of [94] 
and the application of the linear warp evolution model 
in [ 95 ]. We note that while the Bardeen-Petterson effect 
is not manifest in our models, we nevertheless do ob¬ 
serve a global partial realignment of the BH spin with 
the disk angular momentum, caused by the growth of 
the PPL This alignment of the disk angular momentum 
with the BH spin has also been observed in the post¬ 
merger evolution of a tilted accretion disk in [27]. The 
authors conclude that the transport of angular momen¬ 
tum by non-axisymmetric shock waves in the disk might 
be responsible for such a Bardeen-Petterson-like effect. 
Elucidating this issue deserves further studies. 

Finally, to check for Lense-Thirring precession, we plot 
the evolution of the inner region disk twist cr(r) from 
t = 3 orbits onwards for all our tilted models in Fig. [^ 
We see that all models become twisted (as a{r) varies 
with r) in the regions for radii up to 5. The plots further 
indicate solid precession of the disks, as the radial profile 
of cr(r) is increasing smoothly in time almost indepen¬ 
dently of the radius in the outer regions of the disk and 
does not remain 0 for large radii. This is consistent with 
the solid-body precession found in the fixed spacetime 
simulations of [32]. Models D2 actually show a slightly 
growing twist for larger radii. In this trend, the higher 
initial spin models show flatter profiles, while there is 
some radial modulation of the profiles for the lower spin 
simulations. A striking feature of the plots is the much 
smaller difference in accumulated twist for large radii be¬ 
tween different initial spin magnitude for models NCI. In 
models CIB and D2 the higher spin models show a much 
higher twist at large radii than their low spin counter¬ 
parts. As the long survival of the m = I mode in the NCI 
models seems to have an effect on the evolution of the 
twist, we plot this evolution only up to the saturation of 
the PPI for these models in order to properly visualize 
the solid body precession of the disk. 


H. Angular momentum transport 


In this section we investigate the transport of angular 
momentum in the disks. To illustrate the transport of an¬ 
gular momentum during the evolution of our models, we 
plot spacetime t — r diagrams showing the magnitude of 
the angular momentum in radial shells, computed using 
Eq. <§• Such diagrams are displayed in Eig.[^ The mag¬ 
nitude of the angular momentum in each shell, ||J(r)||, 
has been rescaled to the global maximum value attained 


in any shell during the evolution. Eig. 24 shows the evo- 
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FIG. 23: Spacetime t — r diagrams of the radial disk twist profile a{r) for models ClBaOl, ClBaOS (top two panels), NClaOl, 
NClaOS (middle two panels), and D2a01 and D2a05 (bottom two panels). From left to right the panels show the initial tilt angles 
/3o of 5°, 15° and 30°. Each panel shows the radial profile of the inner region disk twist from 3 orbits onwards. 


lution of ||J(r)|| for the three untilted models ClBaObO 
(top panel), NClaObO (middle panel) and D2a01b0 (bot¬ 
tom panel). This figure reveals that compared to the lat¬ 
ter two models, model ClBaObO shows a radical redis¬ 
tribution of angular momentum in the radial direction. 


This drastic difference is due to the development of the 
PPI in model ClBaObO. Following the initial perturba¬ 
tion we see that the contours of ||J(r)|| start growing 
in a wave-like manner, transporting angular momentum 
outwards. This happens up until the saturation of the 
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PPI at about 10 orbital periods. As noted by [96] during 
the growth phase the PPI is thought to be an effective 
mechanism for angular momentum transport in thick ac¬ 
cretion disks, where a right combination of rotation and 
pressure at the boundaries of the disk allows the grow¬ 
ing non-axisymmetric modes to transport angular mo¬ 
mentum outwards. In [84], the authors showed as well 
that the development of nonaxisymmetric instabilities is 
closely connected to the transport of angular momentum. 
The authors observed that disks with a larger q than the 
critical value given in [13], q > = 2 — v/3, were still 

unstable, although the angular momentum transport was 
slower in those cases. After saturation, the new distribu¬ 
tion of angular momentum seems PP-stable and does not 
exhibit drastic changes for the rest of the evolution. 

We can compare the evolution of ||J(r)|| to the de¬ 
velopment of the m = 1 non-axisymmetric mode in the 
disk, which drops to a value of around 1% of its max¬ 
imum value after saturation (see top panel of Fig. 12). 
As described in [96], the development of the PPI cru¬ 
cially relies on the properties of both the inner and outer 
boundaries of the disk. In [97] it was shown that a small 
amount of accretion (which modifies the inner bound¬ 
ary of the disk) was sufficient to saturate the growing 
PPI. Furthermore, wider tori are known to be not as vi¬ 
olently unstable as slender tori [14]. In subsections IV C 
and |IVD| we already showed that the accretion rate and 
the growth of the m = 1 non-axisymmetric mode satu¬ 
rate at the same time, but the t—r diagram of Fig. [^ also 
shows that the outer angular momentum boundary of the 
disk is pushed outwards with growing amplitude as well. 
It therefore seems that both mechanisms are saturating 
the PPI simultaneously in model ClBaOO. 

Model NClaObO, displayed in the middle panel of 
Fig. 24, has a non-constant specific angular momen¬ 
tum distribution, and shows a very different evolution 
of ||T(r)||. The inner region of the disk shows a gentle 
reduction of angular momentum during the entire evo¬ 
lution and no sudden and drastic redistribution is found 
as occurred in model ClBaObO. The outer contour grows 
to encompass the entire outer domain shown in the plot. 
Another difference is the almost complete absence of vari¬ 
ability of the inner region of the disk, which is in contrast 
to model ClBaObO, where oscillatory behavior is seen in 
the inner boundary of the disk too. The specific angular 
momentum profile of model NClaObO is between constant 
and Keplerian, with Keplerian disks being essentially PP- 
stable. Model NClaObO still develops the PPI, albeit in a 
more mild manner than model ClBaObO, with the torus 
of model NClaObO being initially wider. 

Finally, model D2a01b0 (bottom panel), the lightest of 
our m odels, is PP-stable as already shown in subsection 
IV D The corresponding spacetime diagram of ||J(r)|| 


shows no pronounced outward transport of angular mo¬ 
mentum as the evolution proceeds, the profile remaining 
very similar to its initial configuration. The inner region 
with a high angular momentum magnitude is thinned 
out slightly, and the profile spreads a bit overall, but the 
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FIG. 24: Spacetime t — r diagram showing the time evolu¬ 
tion of the radial profile of the angular momentum magni¬ 
tude II J(r)|| for untilted models ClBaObO (top panel), NClaObO 
(middle panel) and D2a01b0 (bottom panel). 




changes are nowhere near as pronounced as for models 
ClBaObO and NClaObO. 

As we have seen in the analysis of the modes, the am¬ 
plitude of the m = 1 mode drops sharply upon the satu¬ 
ration of the PPI in the ClBaOO model. The subsequent 
evolution of the angular momentum profile remains more 
constant, similar to what we observe for the entire evo¬ 
lution of the PP-stable model D2a01b0. There is a clear 
split between the stages of the evolution: the first during 
which the PPI develops and saturates, and the subse¬ 
quent evolution of an essentially PP-stable torus. 

Model NClaOO, on the other hand, shows a very differ¬ 
ent behavior: as already seen in the mode analysis, the 
m = 1 mode amplitude remains roughly at its saturation 
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value for the rest of the evolution, which is lower than 
the values attained by models GIB. The persistent m = 1 
mode is seen to continuously transport angular momen¬ 
tum outwards till the end of the simulation. The restruc¬ 
turing of the disk by the transport of angular momen¬ 
tum is more long-lived, compared to the drastic change 
the saturation of the PPI brings about in model ClBaOO. 
These findings seem to confirm those found in [84]. Us¬ 
ing the same kind of spacetime diagrams we have also 
checked the dependence of the angular momentum trans¬ 
port on the BH spin a and on the initial tilt angle Pq for 
all models of our sample. We find that the angular mo¬ 
mentum transport shows weak dependence with the tilt 
angle (the higher the tilt angle the slowest the trans¬ 
port) and shows essentially no dependence with the spin, 
at least for the moderate values of the BH spin we could 
afford in our study. This is a remarkable result, as it 
demonstrates that models GIB and NGl are PP-unstable 
for a wide range of initial BH spin magnitudes and tilt 
angles. 


I. Gravitational waves 


The recent numerical simulations of m have shown 
that the growth and saturation of the PPI makes BH- 
torus systems strong emitters of large amplitude, quasi- 
periodic gravitational waves, potentially detectable by 
forthcoming ground-based and spacecraft detectors. It 
was found in particular that the m = 1 non-axisymmetric 
structure survives with an appreciable amplitude after 
the saturation of the PPI nonlinear growth, which leads 
to the emission of quasi-periodic GWs with a large am¬ 
plitude. (First estimates through full GRHD simulations 
of the GW detectability by the PPI instability in self- 
gravitating BH-torus systems were obtained in [18] where 
the nonlinear saturation phase was almost reached - see 
also [98] for earlier proposals and alternative estimates 
of the gravitational waves emitted by non-axisymmetric 
instabilities in such systems.) In this section we turn 
our attention to analyze the implications of the PPI in 
our tilted BH-torus systems regarding the emission of 
GWs, namely the dependence of the resulting gravita¬ 
tional waveforms and spectra on the BH spin magnitude 
and initial tilt angle. 


In Fig. 25 we plot the real part of the (/,m) = (2,2) 
mode of the outgoing part of the complex Weyl scalar 
T 4 for models GlBaOl and NGlaOl and for all initial tilt 
angles Pq. This quantity has been computed at an ex¬ 
traction radius r = 640. Fig. [25] shows that, in agreement 
with m, all models display strong emission of gravita¬ 
tional waves and that the emission persists for many dy¬ 
namical timescales well after the saturation of the PPI. 
The peak amplitude for models GlBaOl is about an order 
of magnitude higher than for models NGlaOl. We recall 
that model GIB has a slightly higher initial disk-to-BH 
mass ratio than model NGl (0.16 and 0.11, respectively) 
and that the development of the PPI in models GIB is 
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FIG. 25: Time evolution of the I — m — 2 mode of the real 
part of the Weyl scalar T 4 multiplied by the extraction ra¬ 
dius (r — 640). The top panel shows models GlBaOl and the 
bottom panel models NGlaOl. 
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more pronounced than in the non-constant specific an¬ 
gular momentum models NGl (see Fig. [^. These two fea¬ 
tures are responsible for the difference in the peak ampli¬ 
tudes between the two initial models. We also note that 
the peak amplitude in the GlBaOl models, reached at PPI 
saturation, is much higher than the amplitude of the re¬ 
maining signal, whereas in models NGlaOl the variations 
are not so pronounced. As we showed in Fig. upon 
PPI saturation the dominant m = 1 PPI mode drops 
to about 1% of its peak value in models GIB, while it 
remains at a similar strength for models NGl. 

Fig. [^ also shows that the GW signal has a weak 
dependence with the initial tilt angle, particularly for 
the non-constant specific angular momentum models NGl. 
The smallest peak amplitudes are attained for the most 
tilted BH spacetime (/do = 30°). In any event, the differ¬ 
ences found in the values spanned by the peak amplitudes 
with regard to /3o are not too significant. 

We use the fixed frequency integration (FFI) described 
in [ 88 ] in the integration of the T 4 data to obtain the 
GW strain. In Figs. |26l^ we plot the effective strain 
heff/MPc as a function of frequency for all models, plac¬ 
ing the sources at 1 MPc, in order to see any possible im¬ 
print of the initial spin magnitude and tilt angle on the 
GW spectrum. As described above, we use the same time- 
step in the 7 outermost refinement levels (At = 0.32 and 
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FIG. 26: Spectrum of the effective strain for models GIB, show¬ 
ing the I — 2^m — 1 modes (top) and I — m — 2 modes 
(bottom) for the source located at 1 Mpc. 
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FIG. 27: Spectrum of the effective strain for models NCI, show¬ 
ing the I — 2^nn — 1 modes (top) and I — m — 2 modes 
(bottom) for the source located at 1 Mpc. 


0.64 for higher and lower resolution runs, respectively), 
which means that the time-step at the extraction radius 
of r = 640 is sufficiently small to sample the GW multi¬ 
pole signal with a of 1.28 (2.56) for the higher (lower) 
resolution runs, respectively. Due to the Nyquist crite¬ 
rion, we can therefore resolve the GW spectra up to ~ 80 
(40) kHz. The top panels of Figs. 26 and 27 correspond 
to the (/, m) = (2,1) mode and the bottom panels to the 
(/,m) = (2,2) mode. The highest power is usually found 
at a frequency of about IkHz. For the effective strain of 
the (/, m) = (2,1) mode for models CIB shown in the top 
panel of Fig. we can clearly identify two trends. On 
the one hand there is a difference at high frequencies be¬ 
tween the models with either no or low (a = 0.1) initial 
BH spin (models ClBaOO and ClBaOl) and the high spin 
models (a = 0.3). The spectra for these two groups show 
a different slope from / ^ 2kHz onwards. On the other 
hand, for models ClBaOO and ClBaOl, an increase in the 
initial tilt angle leads to an increase in the power for fre¬ 
quencies / > 2kHz. This also seems to be present at the 
beginning of the spectrum although it is somewhat less 
clear. 


Correspondingly, the power spectra of theJLm) = 
(2, 2) mode shown in the bottom panel of Fig. neatly 
splits the models in three groups, according to the ini¬ 
tial spin magnitude of the BH. We find that the slope 
of the spectra becomes steeper from / ^ 2kHz onwards 


the higher the BH spin. In all cases, the slope of the ef¬ 
fective strain with the frequency shows essentially no de¬ 
pendence with the initial tilt angle for the various groups 
of models characterized by the same BH spin. The power 
in model ClBaOO is at least an order of magnitude larger 
for all frequencies than that of the tilted models. This is 
because the spiraling movement in the equatorial plane 
of both the BH and the m = 1 overdensity PPI “planet” 
emit GWs predominantly in the I = m = 2 mode which 
are emitted in a direction perpendicular to the equato¬ 
rial plane. This is an optimal situation for the model with 
aligned BH and disk spins. However for the tilted mod¬ 
els the “orbit” of the BH and of the overdensity lump in 
the disk is not confined to the equatorial plane anymore 
because of the reaction of the disk to the tilted BH spin. 
This results in a decrease in the GW power. 

The spectra of the non-constant specific angular mo¬ 
mentum models NCI displayed in Fig.is markedly dif¬ 
ferent to those of constant specific angular momentum 
tori. Namely, their dependence on the frequency is sim¬ 
ilar for all modes (same slope) irrespective of the initial 
tilt angle and of the BH spin. No split in two or three 
families is found according to the BH spin. The effective 
strain of the (/,m) = (2,1) mode for models NCI shows 
a plateau in the low frequency part (up to ~ 0.5kHz) for 
models NCla03bl5 and NCla03b30. The spectra shows a 
prominent peak around 0.6kHz for all models but these 
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FIG. 28: Spectrum of the effective strain for models D2, show¬ 
ing the / = 2, m = 1 modes (top panel) and / = m = 2 modes 
(middle panel) and / = 3,?n = 1 modes (bottom panel) for 
the source located at 1 Mpc. 


two. For the I = m = 2 mode we find that the curve 
for the untilted model NClaOO lies more than an order 
of magnitude above the curves corresponding to the ini¬ 
tially tilted disks, as seen also for the ClBaOO model in 
the previous figure. 

Finally, in Fig. we show the effective strains of the 
(/,m) = ( 2 , 1 ) mode (top panel), ( 2 , 2 ) mode (middle 
panel) and (3,1) mode (bottom panel) for models D2. 
In the spectra of the (/,m) = (2,1) mode there is a 
clear variation with the initial tilt and BH spin mag¬ 
nitude, with both increasing the power of the mode for 
all frequencies. Contrary to model CIB and similar to 
model NCI, we find now that model D2 does not show 
the grouping of models depending on the BH spin, and 


the associated frequency dependence (i.e. steeper slopes 
for higher spins). As this model is PP-stable and is built 
with an initially constant distribution of specific angu¬ 
lar momentum, its dynamics is somewhat between that 
of the other two models. The spectra of the I = m = 2 
mode (middle panel) show in particular a very different 
behavior to those shown for models CIB and NCI, namely 
that the model with zero BH spin does not show signif¬ 
icantly larger power than the models with a = 0.1 and 
a = 0.5. Since models D2 do not develop the PPI, there is 
no strong emission of GW in the I = m = 2 mode due to 
the absence of the spiral motion of both BH and overden¬ 
sity ‘planet’ that were observed for the other two initial 
models. Therefore, the spectrum of the untilted model 
D2a01b0 does not show a higher power density than the 
rest of the tilted D2 models. The two models D2a05bl5 
and D2a05b30 both show two peaks in the low frequency 
part of the spectrum, as well as a flatter slope between 
^ 4.5 and ^ 6.5kHz. 

The absence of the PPI in the D2 models results in 
a GW emission where the I = m = 2 mode is not the 
dominant mode. To show this, we plot the effective strain 
of the (/, m) = (3,1) mode in the bottom panel of Fig. 
The power is comparable with that of the (/,m) = (27T) 
and I = m = 2 modes, which is not the case for models 
CIB and NCI that both develop the PPL The plot of the 
effective strain of the (/,m) = (3,1) mode shows that 
same trend as that of the (/,m) = ( 2 , 1 ) mode, i.e. the 
power density increases with spin magnitude and initial 
tilt angle. Furthermore, we can see clear quasi-periodic 
oscillations in the low frequency part of the spectrum 
for models D2a05, as well as a developing plateau which 
becomes fiatter with higher initial tilt angle from ^ 4.5 — 
6.5kHz. 

In order to estimate the detectability of the GWs emit¬ 
ted by these systems, we follow the analysis in [19]. We 
assume that the BH-torus system of models NCI will ra¬ 
diate GWs for the entire accretion timescale tacc (which 
is the time needed for the entire disk to accrete). This 
is because in these models, the m = 1 overdensity lump 
and the BH form a long-lasting quasi-binary system. As¬ 
suming that the accretion rate will remain at the levels 
it attained by the end of the simulations (see Fig. [Tq| ), 
the lowest estimate for the accretion timescale for models 
NCI is tacc ~ 2 X 10^. From this timescale, we can esti¬ 
mate the number of GW cycles at the peak frequency, 
which yields 100 cycles as a lowest estimate for mod¬ 
els NCI. The peak amplitude of the GW strain will then 
be amplified by A^cycies ~ 10 ? while we assume no 
such amplification for neither models CIB (where the PPI 
is rapidly damped after saturation) nor for models D2 
(which are PP-stable). We note that these estimates are 
in very good agreement with the corresponding findings 
reported in m- The corresponding peak amplitudes for 
the BH-torus systems located at 50 Mpc, together with 
the advanced LIGO (aLIGO) sensitivity curve, are plot¬ 
ted in Fig. showing that the GWs emitted by models 
NCI could be detectable due to the long-lasting emission 
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FIG. 29: Plot showing the hypothetical peak amplitudes of 
the effective strain, inferred from the accretion timescale, for 
the three different initial models. The sensitivity curve for 
advanced LIGO with zero-detuning and high-power configu¬ 
ration is shown in black. 



of GWs due to the persistent m = 1 mode. Note that 
the planned factor-of-3 upgrade of aLIGO j99l41Qlj will 
further improve the detectability. 


V. DISCUSSION 


In this paper we have presented the results of an ex¬ 
tended set of numerical relativity simulations of massive 
accretion disks around tilted Kerr black holes. We have 
considered three different thick accretion disks of varying 
mass and specific angular momentum magnitude and dis¬ 
tribution, along with different black hole configurations 
with two spin magnitudes a and four initial tilt angles 
On the one hand the motivation for this work has 
been to extend the investigations of isaisa] in the test- 
fluid approximation to full general relativity, analyzing 
the effects of the self-gravity of the disk on the BH-torus 
dynamics. On the other hand, our work has also served 
to enlarge the parameter space of the existing numeri¬ 
cal relativity simulations HU US] by accounting for tilted 
BH-torus systems for the first time. 

For the models with disk-to-BH mass ratios of 0.044 — 
0.16 we have studied, we have found that the assumption 
of using a fixed background spacetime is not complete as 
we have observed significant precession and nutation of 
the tilted black hole as a result of the disk evolution that 
cannot be accounted for in fixed spacetime simulations. 
For some of our models the precession rate attains a fairly 
constant value by the end of the evolution. The BH nu¬ 
tation rate is also seen to be drastically modulated for 
those models that develop the PPI, showing that the de¬ 
velopment of non-axisymmetric modes in the disk exerts 
a torque on the central BH. 

Our simulations have revealed that the development 
of the PPI seems to be universal for the range of ini¬ 
tial spin magnitudes and tilt angles investigated. Models 


GIB and NCI, both PP-unstable when the central BH is 
non-rotating or untilted, remain PP-unstable when the 
BH is rotating or tilted. The reverse also seems to be 
true: model D 2 , PP-stable for an untilted, non-rotating 
BH remains PP-stable for all spins and initial tilt an¬ 
gles considered. As we mentioned before, model D 2 , hav¬ 
ing the smallest mass ratio in our study, is the most 
likely outcome of BH-NS mergers based on astrophysical 
considerations. Thus, our findings are relevant to gauge 
the practical role of the PPI in thick post-merger ac¬ 
cretion disks around black holes. In agreement with the 
numerical relativity simulations of m the growth of the 
m = 1 PPI mode manifests itself in the formation of 
a counter-rotating overdensity lump (or “planet”) that 
forms a “quasi-binary” with the central BH during its 
existence. This causes the BH to start moving in a spiral 
trajectory for as long as the “planet” exists. For models 
GIB the “planet” disperses quickly upon the saturation of 
the PPI, while in models NCI the overdensity structure 
persists much longer, with the m = 1 mode amplitude 
remaining at its level of saturation. When the PPI devel¬ 
ops around a tilted Kerr BH, the binary motion of BH 
and “planet” is not restricted to the x^-plane but also 
shows some vertical motion. For models GIB in particular, 
the rapid destruction of the non-axisymmetric structure 
causes a mild kick to the BH-torus system in the vertical 
direction which is also imprinted in the time evolution 
of the linear momentum radiated away by gravitational 
waves in the same direction. 

The evolution of the disk around the tilted BH can 
cause a significant twisting (differential precession) and 
warping (spatially varying tilt) in the disk. We have mon¬ 
itored the evolution of the twist cr(r) and tilt z^(r) of ra¬ 
dial shells during the simulations. For models GIB (NCI) 
we have found a phase of rapid (mild) realignment of the 
total angular momentum vector of the disk, Joisk, and 
the BH spin, Jbh, during the growth of the PPI. We at¬ 
tribute this alignment to the development of the m = 1 
structure in the disk. Our simulations have also confirmed 
the presence of significant differential twisting of the disks 
due to Lense-Thirring precession. For all models, the cu¬ 
mulative twist is higher for higher initial BH spins, as 
expected, and the outer regions of the disks precess as a 
solid body. The evolution of the tilt profiles of models D2 
is similar to what [32] observed, with the development 
of a noticeable peak in the innermost region of the disk. 
For models CIB the PPI causes a drastic change in the 
tilt profile upon saturation. Correspondingly, models NCI 
develop strong tilt oscillations in the regions close to the 
BH after the PPI saturates and the m = 1 structure per¬ 
sists till the end of the simulations. By interpreting the 
long-lived m = 1 “planet” as an eccentric distribution of 
matter, we suggest that the so-called Kozai-Lidov mecha¬ 
nism, in which eccentricity and inclination are exchanged 
in an oscillatory manner, might explain the tilt angle os¬ 
cillations we observe, in analogy with the recent findings 
of [92|. In none of our models we observe the Bardeen- 
Petterson effect, as z/(r) 7 ^ 0 in the inner region of the 
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disks. 

The PPI is thought to activate the outward transport 
of angular momentum. This is indeed the case in our 
simulations, as the evolution of the angular momentum 
profiles shows. For model ClBaObO, the angular momen¬ 
tum profile contours start growing in a wave-like manner 
up until the sudden saturation of the PPI. After satu¬ 
ration there is no more transport of angular momentum 
and the evolution resembles that of the PP-stable model 
D2a01b0. For model NClaObO the persistence of the m = 1 
structure in the disk continuously transports angular mo¬ 
mentum outwards for the entire evolution of the disk. 

As [19] showed, the development of the PPI and the 
corresponding overdensity m = 1 lump in the disk cause 
the long-term radiation of gravitational waves, predom¬ 
inantly in the I = m = 2 multipole mode, as expected 
from the radiation emitted by a binary system. Our simu¬ 
lations have confirmed these results, showing that models 
CIB and NCI do indeed radiate mainly in that particular 
mode. The rapid destruction of the m = 1 structure in 
models CIB causes the amplitude of the emitted GW sig¬ 
nal to drop significantly after the PPI saturates, while 
the peak amplitude is closely correlated with the time 
of saturation. Models NCI on the other hand, emit at 
an amplitude similar to that attained at saturation, as 
the m = 1 structure survives until the end of the sim¬ 
ulations. We have also calculated the effective strains of 
the GW signals. While models CIB show different spectra 
for different initial spin magnitudes and tilt angles, there 
seems to be no such trend for models NCI. For those two 
set of models, the GW emission is predominantly in the 
I = m = 2 mode while in models D2 there is also signif¬ 
icant power in the / = 3, m = 1 mode. In addition, the 
strain spectra of models D2 show a clear trend to higher 
values with increasing initial spin and tilt angle. Finally, 
we have briefly touched on the issue of the possible mod¬ 
ulation of the GW signal caused by the BH precession. 
For the fastest precessing models of our sample, models 
ClBaOl, we have shown that we would need about 80 or¬ 
bits for a complete BH precession cycle in order to see 
such an effect, assuming the precession rate remains con¬ 
stant beyond the 20 orbit mark we could afford in our 
simulations. We plan on exploring such a long-term sim¬ 
ulation in the future in order to obtain the complex GW 
signal from the kind of precessing BH-torus systems we 
have explored in this work. 
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Appendix A: The BSSN evolution equations 

The BSSN formulation introduces a set of new auxil¬ 
iary variables for the evolution, defined in terms of the 
ADM quantities: the conformal factor 

0 = In (^det7ij) , (Al) 

the conformal spatial metric 

7ii = , (A2) 

the trace of the extrinsic curvature 

K = (A3) 

the conformal trace-free extrinsic curvature 

, (A4) 

and the conformal connections 

(A5) 

Using these auxiliary variables together with appropriate 
gauge choices leads to a strongly hyperbolic system |lQ2j 
that allows for long-term stable numerical evolutions. 

The evolution equations used in McLachlan are the 
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(A 6 ) 

(A7) 
(AS) 
, (A9) 

(AlO) 

(All) 


(A 12 ) 


(A13) 

where we have used the following shorthand notation, 
do = dt- dj. 

The coupling to the stress-energy tensor is done via the 
usual projections with the spatial metric ^ij and normal 


vector 

E = (A14) 

= T(Too-/?*Toi + /3'yrJ), (A15) 

Sij = (A16) 

S = Si = 7'^' Si^ , (A17) 


Si = -7i^ n, - ^ (^oi - Tij) ■ (AIS) 

This system constitutes the so-called ^-variant of the 
BSSN formulation (see m for other possible varia¬ 
tions). 

McLachlan uses the 1+log slicing for the evolution of 
the lapse [TM]: 


following: 

doa = -a^ - Ko{x^")), 

doK = -YWiDja + aiA^^ Aij + ^K^) 

+ME + f^Sij), 

9o/3* = G{a,4>,x^) , 

doB^ = e-^‘l’H{a,(l),x^^)dor -v\B\a,x^^) 

do4> = -^{aK-dkP>^) , 
dolij = -2 a Aij + 2 ^ % dkp'" , 

d^Aij = (a Rij + aRt - DiD^aj 

+aK Aij —2a Ai)^ A^j + 2 Ay-^idj^P^ 
- 9,nae~^'*‘ , 

doV = -2A^^dja 

+2a{fliA^^ +&A^i K,j 

_fi + ? p 1 f+ 0 

-167ra7®*' Sk , 


f{a,4),xk') = 

a 

Ko{xn = 0 , 


(A19) 

(A 20 ) 


where q{r) is a function that attenuates the T-driver de¬ 
pending on the radius and is a damping parameter. 
We evolve the damping parameter 77 using the formula¬ 
tion of m adapted to a single puncture. This approach 
considers a dynamically evolved damping parameter as 
well as a radial falloff in 77. Assuming a constant gauge 
function 77 causes the shift to diverge at mesh refinement 
boundaries at already very early times. 


Appendix B: GRHD evolution 

In the so-called Valencia formulation that is evolved 
in GRHydro, the evolved variables, called the conserved 
variables^ are defined in terms of the primitive variables 
via the relations 


D = ^pW, 

(Bl) 

Si = ^phW^Vi, 

(B2) 

T = ^{phW^-p)-D, 

(B3) 


where p is the rest-mass density, h is the specific enthalpy, 
h = l + e + p/p, eis the specific internal energy, p is the 
isotropic pressure, and is the three-velocity. In addi¬ 
tion, W is the Lorentz factor and 7 is the determinant of 
the spatial metric. The primitive variables are calculated 
from the conserved quantities at each time step of the 
evolution using a root-finding procedure ESI- 

The corresponding first-order flux-conservative hyper¬ 
bolic evolution system is: 


dU dF^ 
dt dx^ 

with the vector of conserved quantities 
U = [D,Si,T] , 


(B4) 


(B5) 


the vector of fluxes 


F*= [Dv\SjV^ + 5^jP,Tv^ +pv^] , (B 6 ) 


with v* = an* — /3*, and the vector of sources 


0 


S = 


^S’^’^diju + Skdijd'^ - (t -I- D) dia 


aS^^Kki - S^dja 


(B7) 


The stress-energy tensor used in Eq. 0 is that of a 
perfect fluid, given by 


and the following T-driver shift condition ina for the 
evolution of the shift: 


G{a,(p,xk') = 

(A 21 ) 


(A 22 ) 

p\B\a^x^) = q{r) ^ 

(A23) 


= phu^ . (B 8 ) 

The evolution system is closed with an equation of state 
(EOS) which relates the pressure to the other primitive 
quantities. In our simulations the torus is initially de¬ 
scribed by a polytropic EOS 

p = Kp' 


(B9) 
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FIG. 30: Convergence properties: Top panel: dependence on 
grid resolution for the L2-norm of the Hamiltonian constraint. 
Bottom panel: corresponding order of convergence. Model 
ClBaOlbSO has been used for this figure. 

where K is the polytropic constant and F the adiabatic 
exponent. During the evolution we allow for non-isotropic 
changes in the fluid flow such as shocks and for this reason 
we use an ideal gas EOS, 

p=(F-l)pe, (BIO) 

where F is the same adiabatic index in Eq. ( |B9[ ). As in¬ 
dicated in Table [I| we use F = 4/3 for the evolution of 
all our models. A Gamma-law EOS with F = 4/3 is com¬ 
monly adopted to describe a single-component perfect 
gas in the relativistic regime. Our work focuses in par¬ 
ticular on the comparison with previous work on tilted 
accretion disks by [32] , as well as the works of [18] and [19] 
where this value of F was also adopted. For thick disks 
formed by the merger of two NS this EOS may be a poor 
approximation. We are planning on extending our work 
to more realistic EOSs (including the effects of compo¬ 
sition and pressure from trapped neutrinos) such as the 
one described in 

Appendix C: Accuracy and convergence 

This appendix provides details of the accuracy and 
convergence properties of our simulations. In Eig. [30] we 
plot the time evolution of the I/2-norm of the Hamil- 



FIG. 31: Dependence of the irreducible mass of the BH on grid 
resolution for model ClBaOlbSO plotted for the first 4 orbits. 

tonian constraint (top panel) and the convergence or¬ 
der of this quantity (bottom panel) for the first four or¬ 
bits of model ClBaOlbSO. The plot shows the results ob¬ 
tained by using both the canonical (finest grid resolution 
Ax = 0.02) and two lower {Ax = 0.04 and Ax = 0.08) 
resolutions. We follow m\ and define the convergence 
order q as 

, 1 / lll|H|ir-||HE'-|| \ 

log(/) Vlll|H||r'‘-||H||?'''lir 

where / is the refinement factor / = Axiow/Axmed = 
^^med/^^high = 2. Erom Eig.l^wc see that the order of 
convergence shows some variability between a minimum 
of ^ 1 and a maximum of about 3. The spacetime evolu¬ 
tion is fourth order, while the evolution of the hydrody¬ 
namics is only second order accurate and reduces to first 
order in the presence of shocks. We observe higher than 
second order convergence in the I/2-norm of the Hamilto¬ 
nian constraint. This could be ascribed to the spacetime 
evolution which, due to the small disk-to-BH mass ratio 
and to the inconsistent ID setup, could dominate the er¬ 
ror budget in our simulations. Eurthermore, it is difficult 
to assess the order of convergence of our simulations at 
later times once the PPI sets in, since the truncation error 
of the different resolutions could excite the exponential 
growth of the non-axisymmetric modes at different times 
causing the error to diverge. 

The fractional change in the irreducible mass of the BH 
for model ClBaOlbSO is plotted in Eig.[^ The most strik¬ 
ing feature of this figure is the absence of an initial drop of 
the irreducible mass due to the inconsistent choice of ID, 
contrary to the behavior found in the simulations of m- 
The initial drop found by [18] (which is about 8-10%) 
happened in the first orbit and was due to their metric 
blending technique in the ID setup. After that significant 
initial drop, the irreducible mass in m increases again 
to remain at around 97.5% of its initial value during the 
evolution. Such a precise control of the irreducible mass 
conservation may be related to the use of a (quasi) spher¬ 
ical grid. In contrast, in our setup without AH excision. 
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FIG. 32: Evolution of the fractional error in the irreducible 
mass of a vacuum Kerr BH with the same mass as that of 
model CIB and a Kerr parameter of a = 0.3. Results are shown 
for two different resolutions. 

there is no oscillatory behaviour in the evolution of the ir¬ 
reducible mass at the early stages of the evolution, as the 
ID in quasi-isotropic coordinates dynamically adapts to 
the puncture gauge. Nevertheless, the irreducible mass is 
seen to start decreasing in a smooth fashion from O.Storb 
onwards. By the 4 orbital periods mark shown in the fig¬ 
ure, the irreducible mass has decreased to a minimum of 
about 92% of its initial value, to subsequently increase as 
the simulation proceeds as matter and angular momen¬ 
tum accrete onto the BH. Moreover, Fig. [3T] shows that 
this (unphysical) drop in the irreducible mass does not 
converge away with resolution. 

In order to investigate this issue, we have performed 
two additional vacuum simulations of a BH with the same 
mass of that of the CIB models and with an initial Kerr 
parameter of a = 0.3. In Fig. we plot the correspond¬ 
ing fractional error in the irreducible mass for two (finest 
grid) resolutions. Ax = 0.04 and 0.08, for a total time of 
20 orbits. Similar to our results for the BH-torus setup, 
we observe no convergence in the unphysical drop of the 
irreducible mass during the first ~ 10 orbits, while the 
loss of irreducible mass has a smaller slope for the higher 
resolution run in the later stages of the evolution, as ex¬ 
pected. We note, however, that the magnitude of the frac¬ 
tional error in the vacuum tests is much smaller (about 
two orders of magnitude) than that observed in our BH- 
torus runs. 


Unphysical loss of irreducible mass for constraint- 
violating ID has been reported before in the literature 
(see, e.g. |lQ8I411Q] k It has been interpreted as if the ef¬ 
fect of the constraint violations were equivalent to the 
presence of a negative mass in the system. The absorp¬ 
tion of the violations by the BH causes the area of the 
AH, and therefore also the irreducible mass, to decrease. 
As stated in Section [Hlj our ID is manifestly constraint- 
violating and, furthermore, the constraint violations do 
not converge away with resolution. In Fig. we plot the 
initial profile of the Hamiltonian constraint for the vac- 



FIG. 33: Initial radial profile of the Hamiltonian constraint 
for a vacuum Kerr BH and for model ClBaOSbSO along the 
x-axis from the AH outwards. 


uum BH run with Ax = 0.04 and for the corresponding 
BH-torus ClBa03b30 model. We clearly see a large con¬ 
straint violation (compared to the vacuum case) in the 
region where the disk initially resides. By interpreting 
the loss of irreducible mass as caused by the absorption 
of constraint violations, the lack of convergence in the 
irreducible mass drop in our BH-torus systems can be 
explained. As we have also seen in the vacuum tests, this 
transient feature connected to the initial constraint viola¬ 
tions (much smaller in this case) is very long lived. In our 
BH-torus models the loss of irreducible mass only flattens 
out once enough matter has been accreted. Such large 
fractional errors are not reported in simulations where 
the BH-torus system forms after the evolution of a binary 
BH-NS or NS-NS merger, which highlights the undesired 
effect of employing constraint-violating ID. 
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